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Preface 

The  purpose  of  this  study  was  to  examine  the  effects  of 
airframe  vibration  on  the  accuracy  of  an  inertial  navigation 
system.  The  was  done  by  developing  models  for  the  vibration 
and  using  an  existing  system  model  in  a  Monte  Carlo  sim¬ 
ulation  of  the  system  error  equations. 
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Abstract 


This  study  examines  the  effects  of  airframe  vibration  on 
the  accuracy  of  a  strapdown  inertial  navigation  system.  A 
stochastic  model  of  the  system  error  equations  is  included,  as 
are  two  models  of  airframe  vibration.  Software  subroutines 
for  model  implementation  in  SOFE  are  included. 

A  representative  C-130A  flight  profile  was  developed 
using  a  flight  trajectory  generator,  PROFGEN .  The  system 
errors  induced  in  the  inertial  navigation  by  simulating  this 
mission  are  included  as  are  those  caused  by  vibration. 
Vibration  induced  errors  were  found  to  be  very  small  and 
orders  of  magnitudes  smaller  than  those  caused  by  other  error 
sources . 
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STUDY  OF  THE  EFFECTS  OF  VIBRATION  ON  THE 
ACCURACY  OF  AN  INERTIAL  NAVIGATION  SYSTEM 


I .  Introduction 


Background 

Inertial  navigation  systems  have  been  in  use  for  many 
years.  As  the  technology  of  these  systems  has  changed  and 
matured,  many  problems  have  been  solved  and  the  capabilities 
of  inertial  navigation  systems  have  increased  greatly. 
However,  as  the  capabilities  of  inertial  navigation  systems 
improve,  new  missions  emerge  which  put  even  greater  demands 
on  the  system.  Error  sources  which  previously  were  small 
enough  to  be  ignored,  suddenly  become  significant. 

Problem 

The  Military  Airlift  Command  (MAC)  is  tasked  with 
providing  airlift  to  all  services  of  the  United  States  armed 
forces.  In  order  to  meet  the  needs  of  the  services,  MAC 
requires  the  ability  to  fly  long  missions,  including 
low-level  flight,  without  external  aiding  of  their  inertial 
navigat.on  systems.  This  requires  the  inertial  navigation 
system  to  maintain  low  position  and  velocity  errors  which 
they  are  able  to  do  in  the  relatively  benign  environment  of 
high  altitude  cruise.  However,  there  are  indications  that 
the  more  severe  vibrational  environment  at  low  altitude. 
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caused  by  greater  atmospheric  turbulence  and  increased  air 
density,  may  cause  intolerable  degradation  of  inertial 
navigation  system  performance. 

The  problem  undertaken  in  this  study  will  be  to  analyze 
the  effects  of  vibration  upon  the  inertial  navigation  system 
accuracy.  A  representative  transport  mission  profile,  to 
include  both  high  and  low  altitude  flight,  will  be  used  to 
excite  the  low  frequency  system  errors. 

The  probabilistic  approach  undertaken  here  utilizes 
stochastic  models  to  account  for  uncertainit ies  in  the 
inertial  navigation  system.  Stochastic  process  and  modern 
estimation  theories  will  be  used  to  characterize  the  initial 
conditions,  forcing  functions,  and  the  resulting  system 
outputs . 

Scope 

The  focus  of  this  study  will  be  to  identify  the  relative 
effect  of  vibration  on  overall  system  performance.  This 
study  will  use  an  error  model  of  the  Honeywell  SIGN-III 
strapdown  inertial  navigation  system.  The  model  used  is  one 
which  was  developed  by  Widnall  and  Grundy  (Ref  3).  A  flight 
profile  will  be  generated  based  on  the  performance 
characteristics  of  the  Lockheed  C-130A  aircraft.  This 
aircraft  was  chosen  because  it  is  the  most  common  aircraft  in 
the  MAC  fleet.  It  provides  a  more  severe  vibrational 
environment  than  either  the  C-141  or  the  C-5A,  and  some 
flight  test  data  is  available  on  C-130  vibration  levels. 


Development 


The  initial  portion  of  this  study  entails  the 
implementation  of  the  SIGN-III  error  model  to  perform  the 
error  analysis.  The  model  equations  will  be  implemented  in  a 
digital  simulation  program  called  SOFE  which  was  developed  by 
Stanton  H.  Musick  (Ref  S) .  SOFE  will  be  used  to  generate 
error  statistics  via  Monte  Carlo  simulation  of  the  linearized 
first  order  error  equations  of  the  SIGN-III  system.  A  Monte 
Carlo  simulation  was  chosen  instead  of  a  covariance  analysis 
because  it  was  believed  that  the  Monte  Carlo  approach  would 
require  less  computer  resources  (Ref  13)  . 

Three  additional  supporting  computer  programs,  PROFGEN 
(Ref  9),  SOFEPL  (Ref  4),  and  DISSPLA  9,  a  utility  package  of 
subroutines  useful  for  plotting  (REF  13),  will  be  used. 
PROFGEN  is  a  flight  profile  generator  developed  by  Musick  and 
SOFEPL,  developed  by  Musick,  Feldmann,  and  Jensen,  is  a 
postprocessor  for  generating  sample  statistics  and  plots  in 
conjunction  with  SOFE.  DISSPLA  is  a  post  processing/piott ing 
program  which  will  used  to  produce  plots  of  the  SOFEPL 
resul ts . 

While  it  is  difficult  to  ensure  that  the  error  model 
differential  equations  are  correctly  implemented,  some  degree 
of  confidence  in  model  validity  will  be  achieved  by 
determining  the  model  responses  to  various  initial 
conditions.  The  resulting  error  plots  will  be  compared  to 
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those  produced  in  similar  fashion  by  Widnall  and  Grundy  which 
were  published  in  Reference  6. 

A  representative  C-130  flight  profile  will  be  developed 
and  implemented  in  PROFGEN .  The  resulting  aircraft  dynamics 
data  will  be  used  to  drive  the  error  differential  equations 
as  implented  in  SOFE.  AM  error  sources  will  be  turned  on 
except  that  no  vibration  will  be  used.  This  will,  to  some 
degree,  validate  that  the  model  has  been  correctly 
implemented.  Furthermore,  this  information  can  be  used  in  a 
comparison  of  error  magnitudes  which  Cin  help  define  the 
significance  of  the  vibrat icn-induced  errors. 

The  next  step  will  be  to  develop  a  means  of  simulating 
the  vibrational  environment.  Unfortunately,  the  normal 
method  of  using  shaping  filters  driven  by  white  Gaussian 
noise  to  drive  the  error  equations  causes  an  extremely  heavy 
computer  burden.  The  C-130  vibration  environment  includes 
frequencies  as  high  as  459  Hz  and  Shannon's  Sampling  Theorem 
would  require  sampling  at  a  minimum  of  913  Hz.  For  a 
representative  MAC  mission  of  about  8  hours  (28,800  seconds), 
this  would  require  a  minimum  of  25,920,000  samples.  In 
addition,  the  vibration  power  spectral  densities  consist  of 
very  sharp  peaks  at  the  C-130  propeller  blade  passage 
frequency  of  51  Hz  and  its  first  8  harmonics.  To  generate 
this  power  spectral  density  would  require  a  large  bank  of 
shaping  filters,  which  would  require  £  significant  amount  of 
computer  time  per  sample.  The  combination  of  the  two 
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problems  makes  the  normal  method  of  using  shaping  filters  to 
generate  the  vibration  far  too  computationally  burdensome  to 
be  used  to  inject  noise  into  the  system.  To  alleviate  this 
problem,  the  vibration  will  be  simulated  via  two  other 
methods  and  the  results  of  those  methods  will  be  compared. 

In  the  first  method,  each  vibrational  peak  will  be 
represented  by  a  single  sinusoid.  All  vibrational  energy 
which  is  outside  of  these  nine  peaks  will  be  assumed  to 
contribute  little  to  total  system  errors  and  will  be  ignored. 
This  will  result  in  nine  sinusoids  for  each  axis.  The 
resulting  sinusoidal  series  will  be  used  to  drive  the  error 
differential  equations  with  all  other  error  sources  turned 
off,  so  that  the  error  due  solely  to  vibration  can  be 
determined.  A  more  detailed  explanation  of  this  method  will 
be  given  in  Chapter  2. 

The  second  method  of  simulating  the  vibration  will  also 
represent  the  vibration  as  a  series  of  sinusoids,  but  will 
assume  that  sinusoids  of  different  frequencies  will  combine 
so  that  their  cotribution  will  cancel  itself  out  over  time. 

If  this  assumption  is  correct,  the  cross  terms  cross  terms 
can  be  ignored.  As  a  result,  the  number  of  vibration  terms 
can  be  reduced  by  nearly  an  order  of  magnitude.  This  method 
will  b'-  explained  in  greater  detail  in  Chapter  2. 
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II.  Error  Mode  1  Development 


Basic  Error  Differential  Equations 

In  order  to  analyze  the  performance  of  an  inertial 
reference  system  using  linear  estimation  theory,  a  stochastic 
system  error  model  may  be  expressed  in  the  form  of  a  set  of 
linearized  stochastic  first-order  differential  equations. 
These  equations  are  of  the  form: 

x  =  Fx  ♦  Bu  +  w  (1) 


where 


F  =  Fundamental  Matrix 
B  «=  Control  Input  Matrix 
x  *  Error  State  Vector 
u  ■  Deterministic  Forcing  Function 
w  *  White  Gaussian  Driving  Noise 

Britting  demonstrated  that  the  same  set  of  nine  basic 
equations  could  be  used  for  both  strapdown  and  gimbaled 
inertial  navigation  systems  (Ref  1) .  The  set  of  equations 
consists  of  a  nine-by-nine  matrix  system-independent  error 
model  (Ref  6)  augmented  by  system-dependent  error  forcing 
functions.  The  first  three  states  represent  position  errors, 
the  next  three  are  velocity  errors,  and  the  last  three  are 
tilt  errors.  Further  definition  of  these  states  and  their 
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units  are  given  in  Table  1  on  the  following  page.  The  matrix 
for  the  basic  error  model  (Figure  1)  and  the  corresponding 
notation  (Table  2)  are  given  on  the  following  pages.  The 
first  order  differential  equations  for  the  SIGN-III  inertial 
navigation  system  used  in  this  study,  were  developed  by 
Widnall  and  Grundy  (Ref  3). 


Coordinate  Systems 

The  basic  nine  state  error  model  is  usually  implemented 
in  an  east-north-up  navigation  frame.  But,  the  SIGN-III 
sensitive  axes  are  aligned  in  an  aircraft  body  frame  oriented 
down,  nose,  right-wing  (x,y,z).  The  different  coordinate 
systems  makes  it  necessary  to  convert  sensor  noises  derived 
in  the  body  frame  into  equivalent  forces  in  the  navigation 
frame  before  using  them  in  the  system  error  model.  The 
orientation  of  the  two  frames  may  be  seen  in  Figure  2  on  page 
11.  The  coordinate  transformation  matrix  cP  15 : 

D 


where 


r 

ex  Cey  Cez 
r 

nx  Cny  Cnz 
Cux  Cuy  Cuz 


(2) 


Cex‘  sln* 

Cey*  -sin#  cos  <£ 

C  =  cos#  cos  0 

C  it 
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Cnx=  sin  0  cos  <t> 

CnyS  sintf/  sin  6  sin  <t>  +  cos^cosO 
Cnz=  sin  ^  sin  0  cos  tfr  -  sin^cos 
Cuxe  COS^COSd 

Cuy=  cosflsin^sintff  -  sindcostfr 

C  =  cos  0  sin  4>  cos  \p  ♦  sin<£  sintf/' 
=  Roll  Euler  Angle 
6=  Pitch  Euler  Angle 
\J/=  Yaw  Euler  Angle 


The  outputs  of  the  flight  profile  generator,  PROFGEN 
(Ref  9)  are  in  North-West-Up  coordinates.  These  outputs  must 
be  transformed  into  the  body  and  navigation  frames  for  use  in 
the  error  differential  equations.  This  is  done  implicitly  in 
the  simulation  by  equating  the  west  output  of  PROFGEN  to  the 
negative  east  components  in  the  error  model,  i.e.  one  output 
of  PROFGEN  is  west  velocity  which  is  converted  via  the 
following  equation: 


where 

vg  =  east  velocity  in  navigation  frame 

vw  =  west  velocity  in  PROFGEN  (North-West-Up)  coordinates 

Once  specific  force  and  angular  velocity  are  transformed 
into  the  navigation  frame,  they  must  also  be  transformed  into 
the  body  frame  since  the  sensors  are  fixed  in  the  body  frame. 
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Table  1.  State  Definitions  for  the  Fundamental  Error  Matrix 


State 

Meaning 

Units 

x(l) 

Error  in  east  longitude 

radians 

x  ( 2 ) 

Error  in  north  latitude 

radians 

x  ( 3 ) 

Error  in  altitude 

feet 

x  ( 4 ) 

Error  in  east  velocity 

ft/ sec 

x  ( 5 ) 

Error  in  north  velocity 

ft/sec 

X  ( 6 ) 

Error  in  vertical  velocity 

ft/ sec 

x  ( 7 ) 

Error  in  east  attitude 

radians 

X  ( 6 ) 

Error  in  north  attitude 

radians 

X  ( 9 ) 

L. _ 

Error  in  up  attitude 

feet 

c\. 
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p  /cosl.  -  P/ R  cosL  1/R  cosL 


Figure  1.  F  Matrix  of  the  INS  Genera]  Differential  Equations 


Table  2.  Notation  used  in  Figures  1  and  3 


Symbol  and  Value 

Meaning 

L 

Latitude 

Q  *  7.2921151x10-5  rad/sec 

Earth  rotation  rate 

R  *  20925640  ft 

Radius  of  Earth 

g  =  32.12698510 

Magnitude  of  gravity  vector 

< 

< 

< 

Vel  with  respect  to  earth 

e  n  u 

f  .  t  ,  f 
e  n  u 

Components  of  specific  force 

O  -  Q  cos  L 
n 

North  component  of  earth  rate 

Qus  Q  sin  L 

Up  component  of  earth  rate 

—  V  /R 

n 

Components  of  angular 

v  /R 
e 

velocity  of  E-N-Z  frame 

(v  tan  L)  /R 
e 

Vith  respect  to  earth 

w  = 
e 

Pe 

Components  of  angular 

o>  = 
n 

*>xx  +^n 

velocity  of  E-N-Z  frame 

u>  = 
u 

P\i  +  &  u 

with  respect  to  inertial 

kz  = 

v  /R 
u 

F4  2 

«  2  (Q  v  *  Q  \ 
■“n  n 

r  )  +  p 

u  n 

v  /cos^  L 
n 

F43 

*  P  P  +  p 
u  e  rz 

F44 

=  -  k  -  p  tan 
z  e 

L 

F5  2 

=  -2  Q  v  -  p 
n  e  'n 

v  /cos 
e 

2  L 

F5  3 

8  P n  "  k, 

n  u  z  e 

F63 

e  2g/R  -  (pj  ♦ 

kz 

F9  2 

*  u)  *  o  tan  L 
n  Ku 

v's* 


This  is  accomplished  with  a  transformation  matrix  which  is 
the  transpose  of  that  in  Equation  2. 


Altitude  Channel  Mechanization 

The  vertical  channel  of  all  unaided  inertial  reference 
systems  are  inherently  unstable  (Ref  6) .  Typically,  a 
third-order  baro-inertial  loop  is  added  to  stabilze  the 
vertical  channnel,  but  the  SIGN-111  uses  a  second-order 
stabilization  loop  (Ref  6).  The  SIGN-III  altitude  divergence 
control  equation  is  given  in  Equation  (4) : 


u 


vd 


«=  k,  (h 


-  h  r) 
ref 


(h 


-  h  r) 
ref 


(4)  (Ref  6) 


where 


uvd  «  Baro-inertial  aiding  variable 
hc  »  Geocentric  altitude  above  a  reference  sphere 

href  B  Barometric  indicated  altitude 
k^  and  k2  are  the  channel  gains 


The  SlGN-III  uses  the  baro-inertial  aiding  variable  to 
adjust  its  computed  vertical  velocity  and,  thus,  stabilizes 
its  vertical  channel.  However,  an  error  exists  in  this 
mechanization,  since  hc  is  a  geocentric  altitude  and  hre^  is 
a  measure  of  geodetic  altitude.  The  SIGN-III  equations  given 
in  Reference  3  do  not  include  any  conversion  from  geodetic  to 
geocentric  altitude,  and,  as  a  result,  the  computed  geo¬ 
centric  altitude  will  follow  the  geodetic  barometric 
altitude.  The  error  is  modelled,  in  Reference  6,  as  an 
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additional  altitude  reference  error  which  would  be 
implemented  as  a  deterministic  driving  term  in  the  error 
equations  (Ref  6) .  For  the  purposes  of  this  simulation,  it 
will  be  assumed  that  the  implemented  equations  are  corrected 
(or  are  due  only  to  an  error  in  documentation)  by  changing 
the  computations  so  that  geodetic  altitude  is  converted  to 
geocentric  altitude.  Once  this  assumption  is  made,  the  error 
is  eliminated  and  the  following  error  aiding  equations  can  be 
written: 


h 

ref 


h 

ref 


(5) 

(6) 


Equations  (5)  and  (6)  are  inserted  into  the  differential 
equations  for  states  x(3)  and  x(6)  respectively.  The  value  of 
the  gains  and  ^  have  been  selected  so  that  the 

characteristic  equation  has  a  double  pole  at  s  =  -.01  sec, 
which  provides  a  loop  time  constant  of  about  100  sec.  This 
is  the  same  loop  time  constant  as  that  used  in  the  Litton 
Carrier  Aircraft  Inertial  Navigation  System  or  CAINS  inertial 
navigator  (Ref  4).  The  gain  values  are: 


3  x 

10'2 

-1 

sec 

(7 

3  x 

10"4 

-2 

sec 

(8 
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The  S1GN-I I I  uses  three  torque-to-ba lance , 
single-degree-of-f reedom  floated  gas  bearing  Honeywell 
GG-334A9  gyros.  Since  no  test  data  is  available  on  these 
gyros,  Widnall  and  Grundy  based  their  SIGN-III  gyro  model  on 
the  results  of  a  test  of  a  Honeywell  GG-334A16  gyro  (Ref  6). 
This  is  a  low-noise  version  of  the  GG-334  series  gyro  which 
probably  has  better  performance  than  the  gyros  actually  used 
in  the  SIGN-III.  For  the  purposes  of  this  study,  it  will  be 
assumed  that  the  test  results  are  representative  of  GG-334A9 
performance. 

The  gyros  are  mounted  with  their  input  axes  nominally 
orthogonal  and  are  aligned  with  the  platform  x,  y,  2  axes 
(Ref  6) .  For  the  purposes  of  this  study,  these  will  be 
assumed  to  be  nominally  aligned  with  the  aircraft  body  axes. 
This  is  a  simplification  for  analysis  purposes  in  order  to 
reduce  the  computational  burden.  Typically,  a  better  quality 
gyro  is  used  for  the  roll  axis,  or  the  gyros  are  canted 
relative  to  the  roll  axis  in  order  to  distribute  the  roll 
axis  dynamics  (Ref  13) . 

G-insensitive  gyro  drift,  g-sensitive  gyro  drift,  g 
squared-sensitive  gyro  drift,  gyro  scale  factor  errors,  and 
errors  from  gyro  input  axis  misalignment  will  all  be  achieved 
in  the  simulation  by  augmenting  the  basic  nine  state  error 
model  with  additional  states.  The  g-insensitive  gyro  drifts 


will  be  modeled  as  randou.  walk?,  that  is: 


x  *  w  (9) 

The  noise,  w,  will  be  modeled  as  being  white  and  Gaussian 
with  strengths  as  given  on  the  following  page.  The  state 
initial  conditions  will  be  modeled  as  Gaussian  random 
variables.  The  standard  deviations  for  the  initial 
conditions  are  those  used  in  Ref  6  and  are  listed  on  the  next 
page . 

The  remaining  gyro  effects  are  of  the  form: 

x  «=  0  (10) 

That  is,  they  are  modeled  as  random  biases,  and  obtained  as 
the  output  of  integrators  with  zero  input,  and  initial 
conditions  set  as  above  (Ref  6). 

One  additional  simplification  has  been  made  in  modelling 
the  gyro  torquer  scale  factor  errors.  The  GG-334A9  gyros  are 
torque  rebalanced  with  two  torquers  per  gryo.  The  primary 
torquer  is  used  to  maintain  the  gyro  orientaticn,  while  the 
secondary  torquer  is  used  to  compensate  for  known  errors. 

The  primary  torquer  has  a  high  level  mode  in  which  each  pulse 
corresponds  to  2  x  10"12  radians  and  a  low  level  mode  in 
which  each  pulse  corresponds  to  2  x  10-15  radians  (Ref  6) . 

The  torquer  scale  factor  errors  are  likely  to  change 
depending  on  which  mode  is  used  and/or  whether  the  torquing 
is  in  the  positive  or  negative  direction  (Ref  6) .  As  a 
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result,  a  complete  model  of  torquer  scale  factor  errors  must 
include  4  states  per  gyro.  This  study  will  use  only  two 
torquer  scale  factor  error  states  per  gyro  since  the  C-130  is 
a  large  transport  aircraft  and  is  unlikely  to  change  its 
attitude  rapidly.  As  a  result,  the  high  level  torquing  mode 
will  be  used  rarely,  if  at  all. 

Accelerometer  Error  Model 

The  SIGN-III  uses  three  single  axis  rebalancing 
pendulous  accelerometers.  They  are  mounted  in  a  nominally 
orthogonal  configuration  along  the  x,  y,  and  z  axis  of  the 
platform  (Ref  6) .  The  accelerometer  biases  are  modeled  as 
random  walks.  As  was  done  in  Reference  6,  the  state  :nitial 
conditions  will  be  set  as  Gaussian  random  variables  with 
standard  deviations  as  given  in  Table  4,  as  are  the  strengths 
of  the  white  Gaussian  driving  noises. 

The  SIGN-III  error  model,  from  Reference  6,  also 

augments  the  basic  nine-by-nine  matrix  with  additional  states 
for  accelerometer  scale  factor  error  and  accelerometer  input 

axis  misalignment.  These  are  modeled  as  random  constants. 

The  initial  conditions  are  modeled  as  zero-mean  Gaussian 

variables,  the  standard  deviations  of  which  are  given  on  the 

preceding  page. 

Gravity  Error  Mode  1 

The  SIGN-III  gravity  model  equations  are: 
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4 


(12) 


Gn  =  3uJ2  sin  Lc  cos  Lc/rc 
Ge  =  0  (13) 


where 


J 


2 


L 


c 


eDown,  north,  east  components  of  gravity 
=  Geocentric  radius 
=  Gravitational  constant 
*  3.986005  x  1014m3/sec2 

=  Oblateness  coefficient  of  gravitational 
expansion 

=  Geocentric  latitude 


The  north  equation  has  an  error  in  it.  A  correct  expression 
is : 

Gn  =  -3uJ2R2sin  Lccos  Lc/i£  (13) 

where  R  is  the  equatorial  earth  radius. 

The  SIGN-I I I  equation  has  the  wrong  sign  and  the  omission 
2 

of  the  R  term  causes  the  expression  to  have  an  extremely 
small  value  (Ref  6) .  Probably,  this  is  an  error  in  the 
documentation  rather  than  an  error  in  the  equations  actually 
used.  If  these  errors  are  present  in  the  equations  used, 
they  would  cause  an  error  in  north  velocity  which  could  be 
modeled  as  a  deterministic  forcing  function  as  follows: 


u 


5 


0.00162  g0  sin  2I,C 


(15) 


=  Non-zero  forcing  function  element 
10 


where 


go  =  Equatorial  magnitude  of  gravity  for 
the  reference  ellipsoid 
Lc  *  Geocentric  latitude  (Ref  6) 

This  was  done  in  Widnall  .end  Grundy,  but  in  this  simulation, 
it  will  be  assumed  that  the  correct  mechanization  is  used  and 
no  determinsitic  driving  terms  are  required. 

The  error  model  also  includes  three  states  to  model  the 
effects  of  gravity  deflections  and  anomaly.  All  three  states 
are  modeled  as  first  order  Gauss-Markov  processes  with  the 
correlation  time  being  derived  from  the  aircraft  velocity  and 
the  correlation  distance.  This  is  shown  in  the  equation: 


-(v/d)  x  *  w 


where 


(16) 


x  «=  Error  state 
v  *  Aircraft  velocity 
d  *  Correlation  distance 
w  =  White  noise  of  strength  Q 
Q  =  2o2v/d 

Reference  6  gives  the  gravity  variation  model  for  the  western 
half  of  the  United  States.  These  figures  are  shown  in  the 
following  table: 
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Table  3.  Model  for  Gravity  Variations  (Ref  3) 

Gravity  Error 

- 1 _ 

Standard  Deviation 

Distance 

East-West  Deflection 

26  g 

10  nm 

North-South  Deflection 

17  g 

10  nm 

Anomaly 

35  g 

60  nm 

Barometric  Altimeter  Error  Mode  1 

This  simulation  will  model  two  sources  of  error  in  the 
barometric  altimeter.  These  errors  are  scale  factor  error 
due  to  non-standard  temperature  and  the  error  due  to 
variation  in  altitude  of  a  constant  pressure  surface.  The 
non-standard  temperature  error  will  be  modeled  as  a  random 
constant  scale  factor  effect  with  a  0.03  standard  deviation. 
The  variation  of  a  constant  pressure  surface  will  be  modeled 
as  a  first  order  Gauss-Markov  process.  The  correlation 
distance  for  the  process  is  250  nm  and  the  standard  deviation 
is  500  ft  (Ref  6 ) . 

Truth  Model  for  Zero  Vibration 

The  complete  truth  model,  as  developed  by  Widnall  and 
Grundy  in  Reference  6,  is  a  50-state  model  of  the  form: 

x=Fx+u+w  (17) 

The  states,  x,  are  defined  in  Table  5.  The  fundamental 
matrix,  F,  is  shown  in  Figure  3  with  its  entries  being 
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Figure  3c.  Submatrix  F2 
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Figure  3d.  Submatrix  F3 
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Table 

5.  Error  States  and  Initial  Standard  Deviations 

State 

Definition 

a 

.  x(l) 

Error  in  east  longitude 

-2 

5.7735  x  10  arc  min 

x  ( 2 ) 

Error  in  north  latitude 

—2  1 

5.7735  x  10  arc  min  1 

x  ( 3 ) 

Error  in  altitude 

11.45  ft 

x  ( 4 ) 

Error  in  east  velocity 

1  ft/sec 

x  ( 5) 

Error  in  north  velocity 

1  ft/sec 

x  (6) 

Error  in  vertical  velocity 

0.1  ft/sec 

x  ( 7 ) 

Error  in  east  attitude 

Equation  (36) 

x  ( 8 ) 

Error  in  north  attitude 

Equation  (38) 

x  <  9  > 

Error  in  up  attitude 

Equation  (40) 

X(10) 

X  gyro  drift  rate 

0.025  deg/hr 

X(ll) 

Y  gyro  drift  rate 

0.180  deg/hr 

x  (12) 

Z  gyro  drift  rate 

0.180  deg/hr 

x  (13) 

X  gyro  input  axis  g-sensitivity 

0.2  deg/hr/g 

x  ( 14 ) 

X  gyro  spin  axis  g-sensitivity 

0.2  deg/hr/g 

x  (15) 

Y  gyro  input  axis  g-sensitivity 

0.2  deg/hr/g 

x  (16) 

Y  gyro  spin  axis  g-sensitivity 

0.2  deg/hr/g 

x  ( 1 7 ) 

Z  gyro  input  axis  g-sensitivity 

0.2  deg/hr/g 

x  ( 18 ) 

Z  gyro  spin  axis  g-sensitivity 

0.2  deg/hr/g 

x  ( 19 ) 

X  gyro  input  axis  g  -sensitivity 

0.07  deg/hr/g 

X  ( 20 ) 

Y  gyro  spin  axis  g  -sensitivity 

i 

0.07  deg/hr/g 

X  ( 2 1 ) 

!  Z  gyro  input  axis  g  -sensitivity 

0.07  deg/hr/g 

x  ( 22) 

i 

1 

;  *  9yro  pos.  scale  factor  error 

70  ppm 

Table  5.  (Continued) 


State  !  Definition 


x  (23) 

X  gyro  neg.  scale 

factor  error 

x (24)  ! 

Y  gyro  pos.  scale 

factor  error 

j 

X  (  2  5  )  | 

Y  gyro  neg.  scale 

factor  error 

X  (26 ) 

Z  gyro  pos.  scale 

factor  error 

x  (27)  : 

!  Z  gyro  neg.  scale 

factor  error 

x  ( 28 ) 

i 

i  X  gyro  input  axis 

1  about  Y 

misalignment 

x  (29) 

X  gyro  input  axis 
about  Z 

misalignment 

x  (30) 

Y  gyro  input  axis 
j  about  X 

misalignment 

X  (31) 

i  Y  gyro  input  axis 
.  about  Z 

misalignment 

x  ( 32 ) 

;  Z  gyro  input  axis 
about  X 

misalignment 

x  (33) 

Z  gyro  input  axis 
about  Y 

misalignment 

x  ( 34 ) 

i 

i  X  accelerometer  bias 

x  ( 35 ) 

1  Y  accelerometer  bias 

j 

x  ( 36 ) 

;  Z  accelerometer  bias 

x(37) 

X  accelerometer  scale 

factor  error 


29 


Table 

5.  (Continued) 

State 

• 

Definition 

a 

-  —  ■ 

*(38) 

Y  accelerometer  scale 
factor  error 

3  5  ppm 

x  ( 39 ) 

Z  accelerometer  scale 
factor  error 

35  ppm 

X  ( 40) 

X  accelerometer  input  axis 
misalignment  about  Y 

10  arc  sec 

X  (41) 

X  accelerometer  input  axis 
misalignment  about  Z 

10  arc  sec 

X  (42) 

Y  accelerometer  input  axis 
misalignment  about  X 

10  arc  sec 

x  (43) 

Y  accelerometer  input  axis 
misalignment  about  Z 

10  arc  sec 

X  (44) 

2  accelerometer  input  axis 
misalignment  about  X 

10  arc  sec 

X  ( 45) 

Z  accelerometer  input  axis 
misalignment  about  Y 

10  arc  sec 

X  (  4  6 ) 

Barometric  error  due  to 
variation  in  altitude  of  a 
constant  pressure  surface 

500  ft 

x  ( 47) 

Barometric  scale  factor  error 

0.03 

x  ( 48) 

East  deflection  of  gravity 

26  nq 

x  (49) 

North  deflection  of  gravity 

17  nq 

x  ( 50) 

Gravity  anomaly 

35  tiq 
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Table  6.  SIGN-III  Noise  Density  Matrix,  Non-2ero  Elements 


(Ref  6)  j 


Diagonal 

State 

Noise  Density 

Element 

Variable 

6 

Error  in  vertical  Velocity 

0.045  v/250  nm 

10 

X  Gyro  Drift  Rate 

2 

(0.03  deg/hr)  /hr 

11 

Y  Gyro  Drift  Rate 

2 

(0.02  deg/hr)  /hr 

12 

Z  Gyro  Drift  Rate 

2 

(0.02  deg/hr)  /hr 

34 

X  Accelerometer  Bias 

(10  g) 2/hr 

35 

Y  Accelerometer  Bias 

2 

(10  g)  /hr 

36 

Z  Accelerometer  Bias 

(10  g) 2/hr 

46 

Barometric  Pressure 

2  (500)  v  / 2 5 0  nm 

48 

East  Gravity  Deflection 

2 

2  (26  g)  v/10  nm 

49 

North  Gravity  Deflection 

2 

2  (17  g)  v/10  nm 

50 

Gravity  Anomily 

2 

2  (35  g)  v/60  nm 

Off-Diagonal  State  Noise  Density 

Element  Variable 

N,  =  N  Vertical  Velocity  150  v/250  nm 

6,46  46,6 


Table  7 


SIGN-111  Error  Source  Statistics 


1 _  _ 

State  Variable 

Noise  Spectral  Density 

1 

Random  Walks 

10 

(0.03  deg/hr)2/hr 

11,12 

1 

(0.02  deg/hr)2  /hr 

|  34,35,36 

1 

( 10j/g}  2  /hr 

j - - - - - - - - 

!  First  Order  Markov  Processes 

1 

1 —  _ 

Correlation  Distance 

46 

250  n.m. 

48 

10  n.m. 

49 

10  n.m. 

|  50 

60  n.m. 

! 
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explained  in  Tables  2  and  5.  The  u  vector,  which  is  composed 
of  the  deterministic  driving  terms,  has  only  two  non-zero 
elements  which  are  due  to  mechanization  errors.  In  this 
simulation,  it  is  being  assumed  that  these  errors  have  been 
corrected  and,  as  a  result,  all  driving  terms  will  be  set  to 
zero.  The  strengths  of  the  white  Gaussian  driving  terms,  w, 
are  given  in  Table  6.  This  model  will  be  used  as  the 
reference,  or  control,  against  which  a  similar  model,  modified 
by  the  addition  of  airframe  vibration,  will  be  compared. 

Vibration  Environment 

In  April  of  1982,  the  4950th  Test  Wing,  located  at 
Wright-Patterson  AFB ,  Ohio,  conducted  a  vibration  flight  test 
of  palletized  modular  electronics  rack.  The  purpose  of  the 
test  was  to  establish  actual  vibration  levels  at  various 
locations  on  the  pallet  during  C-130A  ground  and  flight 
operations  (Ref  12).  While  it  is  true  that  this  pallet  is  not 
identical  to  those  used  with  operational  palletized  inertial 
navigation  systems,  it  was  designed  for  use  in  the  0130 
aircraft  and  meets  military  specifications  for  use  with 
airbotne  electronics  equipment  (Ref  12).  For  the  purposes  of 
this  study,  it  will  be  assumed  that  the  test  data  contained  in 
Reference  11  is  accurate  and  representative  of  that  to  which 
an  inertial  navigation  system  would  be  subjected  when  used  on 
board  a  C-130A  aircraft. 

The  vibration  data  used  in  this  study  are  shown  in 
Figures  4,  6,  and  8.  The  energy  levels  are  shown  as  power 
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Figure  6.  Y-Axis  Measured  Vibration.  (Ref  12) 
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Figure  7.  Y-Axis  Model  Vibration 
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Figure  8.  2-Axis  Measured  Vibration  (Ref  12) 
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Figure  9.  Z-Axis  Model  Vibration 
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spectral  densities  with  units  of  g-squared  per  Hz.  As  can  be 
seen  in  the  figures,  the  majority  of  vibrational  energy 
occurs  at  the  fundamental  and  the  first  eight  harmonics  of 
the  51  Hz  blade  passage  frequency.  In  order  to  simplify  the 
vibration  model,  it  will  be  assumed  that  the  vibration 
outside  of  these  nine  peaks  contributes  little  to  system 
errors  and  will  be  ignored. 

One  limitation  of  the  available  data  is  that  no 
information  is  available  about  the  degree  of  correlation 
between  the  vibration  levels  in  each  axis.  Although  it  seems 
likely  that  the  vibration  is  correlated  to  at  least  some 
degree,  in  the  absense  of  any  empirical  data  to  the  contrary, 
it  will  be  assumed  that  the  vibration  in  each  axis  is 
uncorrelated  with  that  in  the  other  two  axes.  This  will 
probably  result  in  some  error  in  the  simulation,  but  it  is 
not  anticipated  that  the  error  will  be  excessive. 

A  good  model  for  vibration  is  provided  by  a  second 
order  Gauss-Markov  process:  the  output  of  a  second-order 
system  driven  by  white-Gaussian  noise  (Ref  8)  The  general 
form  of  the  of  the  output  power  spectral  density  of  this 
filter  is  as  follows: 

$  =  (a2w2*  b2)/(<A  2  w2(2w2-1)  ^2+  w4)  (18) 

This  output  can  be  generated  by  passing  a  stationary  white 
Gaussian  noise  of  strength  Q=1  through  a  second  order  system 


with  the  following  transfer  function: 

G (s)  =  (as  +  b)/(s2  +  2uns  +w2) 

The  system  state  equations  are  given  below 
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in  matrix  form. 

w(t)  (20) 


The  values  for  a,  b/  c,  u>n»  an<^  f  are  selected  to  fit  the 
empirical  data  (Ref  8)  . 

Unfortunately,  this  method  of  modeling  the  vibration 
requires  two  additional  states  for  each  of  the  nine  peaks  in 
each  of  the  three  axis.  Thus,  to  model  the  vibration  via 
second  order  Gauss-Markov  shaping  filters  would  add  54 
additional  states  to  an  already  large  50-state  system  model. 
The  sampling  rate  must  also  satisfy  the  Shannon  Sampling 
Theorem,  which  would  require  that  the  sampling  rate  be  at 
least  918  Hz.  As  a  result,  this  model  of  the  vibration  is 
not  acceptable  for  use  in  this  study. 

In  order  to  reduce  the  computational  burden  caused  by 
the  high  frequency  content  of  the  vibration,  it  was  necessary 
to  develop  a  less  computationally  demanding  method  of 
modelling  the  vibration.  Since  the  vast  majority  of  the 
energy  in  the  vibration  environment  is  found  in  nine  peaks, 
located  at  the  fundamental  and  the  first  eight  harmonics  of 
the  blade  passage  frequency,  it  was  assumed  that  the 


38 


vibration  was  caused  by  the  rotation  of  the  propellers. 

Under  these  conditions,  a  reasonable  model  generates  the 
vibration  as  a  series  of  sinusoids,  one  sine  wave  for  each  of 
the  nine  spectral  peaks  for  each  axis  (Ref  2) .  If  each  of 
the  nine  sinusoids  is  then  treated  as  being  independent  of 
the  other  eight,  the  vibration  for  each  axis  can  be  expressed 
as : 


V(t)  =  A  COS  (ut+  (  ) 


(Ref  2,  Ref  5)  (21) 


where 

V  *  the  vibration  for  a  given  axis 
A  =  square  root  of  twice  the  magnitude  of  the 
corresponding  spectral  peak 
w  *  the  blade  passage  frequency 
<  «  a  randomly  generated  phase  angle. 

Figures  4  through  9,  which  are  on  the  following  pages,  show 
the  power  spectral  densities  of  the  vibration  levels  measured 
in  Ref  12  and  those  generated  by  the  method  described  above. 
It  can  be  seen  that  the  model  is  not  a  perfect  representation 
of  the  actual  environment.  However,  the  majority  of  the 
energy  is  contained  in  the  nine  spectral  peaks  and  there  is 
no  discernible  difference  between  the  peaks  produced  in  the 
actual  environment  and  those  produced  by  the  model  in  either 
frequency  or  magnitude. 
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Since  it  is  unlikely  that  the  phase  angles  will  remain 
constant  for  the  entire  flight,  the  phase  angle  will  be 
shifted  by  a  randomly  generated  phase  angle  every  20  samples. 
The  angle  will  be  generated  by  the  GAUSS  subroutine  in  SOFE 
and  a  standard  deviation  of  one  degree  will  be  used  (Ref  11). 
Both  the  rate  and  standard  deviation  of  the  phase  shifts  were 
set  arbitarily  due  to  a  complete  absense  of  data  on  which  to 
base  them. 

The  computational  burden  can  be  further  reduced  by 
splitting  the  model  into  two  parts.  One  part  will  include 
the  full  fifty-state  error  model  with  all  driving  noises 
except  vibration.  This  part  can  be  propagated  in  SOFE  with  a 
relatively  long  sampling  period  since  its  associated  driving 
terms  do  not  include  any  high  frequency  noises.  The  second 
part  will  include  only  the  vibration  terms  and  those  states 
influenced  by  the  vibration,  i.e.  the  9  states  of  the  Pinson 
model  and  states  13  through  21,  a  total  of  18  states.  The 
vibration  enters  the  model  through  the  g-sensitive  gyro  drift 
states  (states  13  through  18)  and  the  g-squared  sensitive 
gyro  drift  states  (states  19  through  21).  These  states,  in 
turn,  affect  the  attitude  error  states;  7,  8,  and  9.  The 
attitude  errors  are  propagated  through  the  normal  systems 
dynamics,  and,  subsequently,  affect  the  velocity  and  position 
error  states. 

Due  to  the  presence  of  high  frequency  noise  in  the 
18-state  model,  a  short  sampling  period  is  required.  On  the 
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other  hand,  the  50-state  model  no  longer  includes  any  high 
frequency  terms,  and  the  sampling  rate  may  be  reduced  by 
several  orders  of  magnitude.  A  similar  study  used  a  sampling 
period  of  2  seconds  (Ref  13) ,  and  that  was  used  in  this  case. 
While  reducing  the  required  sampling  rate  from  nearly  1000  Hz 
to  0.5  Hz  (a  reduction  by  a  factor  of  nearly  2000)  will  not 
result  in  a  proportionate  drop  in  CPU  time,  nonetheless,  the 
time  required  will  drop  considerably.  Since  the  50-state 
simulation  does  not  need  to  be  repeated  everytime  a  new 
vibration  model  is  simulated,  the  CPU  time  required  is 
further  reduced.  These  savings  more  than  compensate  for 
having  to  simulate  2  separate  models. 

A  further  reduction  in  the  computational  burden  can  be 
made  by  noting  that  vibration  is  zero-mean  and  symmetrical, 
and  the  g-sensitive  gyro  drift  coefficients  are  constant.  As 
a  result,  the  combination  of  the  vibration  and  the  g-sen¬ 
sitive  gyro  drift  terms  will  have  minimal  effect  on  the 
system  attitude  errors.  For  the  purposes  of  this  study,  it 
will  be  assumed  that  eliminating  the  g-sensitive  terms  from 
the  vibration  model  will  result  in  negligible  errors  in  the 
results,  and  they  will  not  be  included  in  the  vibration 
model.  The  resulting  vibration  model  consists  of  the  basic 
nine-state  error  model  with  the  following  terms  driving  the 
derivatives  of  states  7,  8,  and  9: 

V  -cexfyfzDx  'ceyfxfzDl'  'cezfxfyDz 
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During  the  remainder  of  the  study,  this  method  will  be 
referred  to  as  the  sinusoidal  series  vibration  model.  Since 
both  parts  of  the  system  error  models  are  linear,  the 
principle  of  linear  superposition  allows  the  following 
equations  to  be  used  to  find  the  error  state  means  and 
covariances  for  the  system  as  a  whole: 

E  (X)  =  E  (XI)  ♦  E  (X2 )  (25) 

E  (X2)  =  (E  (XI2)  +  E  (X22) )  (26) 

This  method  of  modeling  the  vibration  is  more  efficient 
than  using  shaping  filters,  but  it  still  is  not  efficient 
enough  since  it  requires  performing  51  additions,  75 
multiplications,  and  evaluating  27  trigonometric  functions  in 
order  to  inject  the  vibration  for  each  sample  and  the 
sampling  rate  must  be  fast  enough  to  satisfy  Shannon's 
sampling  theorem,  or  a  minimum  of  918  samples  per  second.  A 
more  efficient  model  is  required. 

Experimental  Integrated  Effect  Vibrat ion  Model 

The  majority  of  the  SIGN-III  error  model  is  made  up  of 
slowly  varing  states  with  only  the  vibration  changing  at  a 
rapid  rate.  For  short  periods  of  time,  the  state  transition 
matrix  can  be  treated  as  being  time  invariant.  If,  during 
this  short  period  of  time,  the  cumulative  effects  of  the 
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vibration  can  be  determined  and  simulated,  then  it  is  not 
necessary  to  generate  the  vibration;  instead  the  effects  can 
be  injected  into  the  rest  of  the  system  at  periodic 
intervals . 

If  the  assumptions  made  in  modeling  the  vibration  as  a 
series  of  sinusoids  are  valid;  then  the  following  equation 
expresses  the  relationship  between  the  east  attitude  error 
state  and  the  vibration: 


vn.x  _  _  r  dx  £  /a  cos  (wit  +a.  )  )  £  (B-  cos  +fi. )) 

X(t)  Lex  i=l  '  i  i  3*1  3  D 

9  9 

r  ,  ,  .  >  r  /n  cos  (wjt  ♦ 8  ) ) 

+  r  Dy  **  (p  cos  wkt  *y  ))  ~ 

Ley  *  k=l  (^k  k  3*1  J  J 

9  9 

r  /  i  ...  F  (A  cos  (wit  *a  )  )  (27) 

-  CezDz  kti  <Ck  COS  <"kt  +V,i5l  1 


where 

X  =  east  attitude  error 

Dx,  Dx,  and  Dx  =  g-square  sensitive  gyro 
drift  coefficients 

A,  B,  and  C  =  the  magnitude  of  the  sinusoids 
representing  the  vibration  in  the  x,  y,  and 
z  axes,  respectively. 

If,  to  use  as  an  example,  the  equation  is  reduced  to  show 
only  the  effects  of  vibration  in  two  axes,  on  the  third  the 


equation  becomes: 
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If  this  equation  is  further  reduced,  again  for  use  as  an 
example,  to  include  only  two  of  the  spectral  peaks  at 
arbitrary  frequencies,  and^,  becomes: 

X(t)  «  -  c  Dx  (A  cos  (w  t  *<*.>  ♦  A,  cos  (w  t  +a  ) ) 

0  X  X  XX  £  4 

(B  COS  (  U!  t  +  3  )  +  B  cos  (a)  t  +  R  )  )  (29) 

1  1  1  P1  2  2  h2 

where 

w  and  u;^  =  integer  multiples  of  the  propel  lor  blade 

passage  frequency 

This  can  be  rewritten  as: 

x  ( t )  *  -  c  Dx  (A  (cosoj  t  cosa  -  sinu  t  sina  ) 
e  x  1  1  1  1  1 

+  p.  (cosco  t  cosa  -  sinu  t  sinoe  ) ) 

2  2  2  2  2 

(cos t  cos/3^  -  sinu^t  sin^) 

+  b2  (cosu^t  cos^  -  sinu^t  sin^))  (30) 


After  the  terms  are  multiplied,  this  becomes: 

2 

X  ( t )  *  -c  Dx  (A  B  (cos  u>  t  cosa  cos Q  -  sin  t  cos  t 
'  ex  11  1  1  HL  1  1 
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(sina1  cos^1  +  cosa^  sin^)  ♦  sin  u^t  sina^  sin/5^  ) 


+  a  B  (cosu t  cosu t  cosa  cosH  -  sinw  t  cosu t 
12  1  2  1  h2  1  2 


sinc^  cosff  -  sinu^t  cosu^t  cosa^  sin^  +  sinu^t 
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sinui  t  sine*  sinrt  )  ♦  A  B  (cosu>  t  cosu  t  cosc* 
sin  2  1  p2  2  1  2  1  2 
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cosa^  sinjg^  +  sint^t  sinai^t  sin^  sin^  ) 

2 

♦  A2  B2  ^C0S  ^2^  cosa2  cos^2 

-  sina^t  cosu^t  (sina2  cos^  +  cos<*2  sin^  ) 


+  sin  u^t  sino<2  sin^  )) 


This  expression  can  be  solved  to  find  X(T)  b^  integrating 
both  sides  of  the  equation.  If  the  limits  of  integration  are 
selected  so  that  the  time  interval  is  an  integer  multiple  of 
pi  divided  by  51  (the  blade  passage  frequency),  X(T)  is  as 
fol lows : 

x (T )  *  (-C  Dx  T/2)  (A  B  (cosa  cos fl  +  sina  sinfl  ) 

6  X  XI  XX  XX 

+  A2  B2  (cosa2  cos/?2  +  sina2  sin^2  ))  (32) 


where  T  ®  the  time  interval 

By  using  trigonometric  identities,  this  can  be  expressed  as: 


X(T)  =  (-C  Dx  T/2)  (A  B  (cos  (a  -B  )  +  cos  (a  +5  ) 
1  '  '  ^ex  11  1*1  11 

♦  cos  {ct^~f3^  ~  cos  (a 

♦  A2  B2  (cos  (<*-£)  +  cos  (a2+^2)  +  COS  (a2~^2) 
-  cos  (C <2  *f^2  ^  *  * 
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This,  in  turn,  reduces  to: 

X(T)  =  “CexDx  T  <AX  Bj_  cos  ♦  A2  B2  cos  (a2~02))  (34) 

Thus,  it  can  be  seen  that  if  T  is  selected  small  enough 
that  the  state  transition  matrix  can  be  assumed  constant  over 
the  time  interval,  and  if  T  is  selected  to  be  an  integer 
multiple  of  pi  divided  by  the  blade  passage  frequency,  then 
this  method  of  simulating  the  effects  of  the  vibration  should 
produce  results  of  acceptable  accuracy.  In  order  to  check 
the  accuracy  of  this  method,  short  Monte  Carlo  simulations  of 
each  method  will  be  accomplished.  The  results  of  these  runs 
will  be  compared  in  Chapter  IV. 

Initial  Conditions 

The  initial  conditions  for  the  first  nine  states  will 
be  based  upon  an  assumed  ground  alignment  at  a  random 
heading.  It  will  also  be  assumed  that  the  baro-inertial 
vertical  channel  has  reached  a  steady  state  condition. 

The  initial  longitude  and  latitude  errors,  *o(l)  and 
xq(2),  will  be  based  on  the  assumption  that  the  navigator  can 

enter  the  alignment  position  to  within  an  error  of  +  /-  0.1 

arcminute.  Since  the  variance  of  a  variable,  distributed 

uniformly  over  a  range  T,  is  T  / 1 2 ,  the  square  root  of  the 

variance  (or  standard  deviation)  may  be  calculated  with  the 

following  equation: 

d=  (0.2  arc  min)  A/T2  =  5.7735  x  10~^  arc  min  (Ref  13)  (35) 
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During  the  Monte  Carlo  simulation,  the  initial  conditions  of 
x(l)  and  x(2)  will  be  generated  by  a  zero-mean  Gaussian 
distributed  variable  based  on  the  standard  deviation  given  in 
Equation  (34).  A  Gaussian  distribution  will  be  used  instead 
of  a  uniform  distribution  since  a  means  of  generating  a 
Gaussian  sample  is  already  available  in  SOFE  and  the 
resulting  error  should  be  minimal. 

The  initial  condition  for  the  altitude  error,  x(3)  , 
will  be  generated  based  upon  the  assumption  that  the  aircrew 
will  be  able  to  enter  the  altitude  to  within  20  feet.  This 
will  be  approximated  by  a  Gaussian  distributed  variable  with 
a  standard  deviation  of  2(20)/vT2  or  11.45  feet  (Ref  13). 

The  initial  standard  deviations  of  the  velocity  states, 
x(4),  x(5),  and  x(6),  will  be  those  recommended  in  Reference 
6.  These  are: 

cr(4)  *  1.0  ft/sec 
a  (5)  =  1.0  ft/sec 
a  (6 )  *=  0.1  ft/ sec 

The  initial  condition  for  the  vertical  velocity,  X  ( 6 )  ,  is 
more  accurate  than  the  other  velocity  initial  conditions  due 
to  barometric  altimeter  aiding. 

The  east  and  north  attitude  errors,  x(7)  and  x(8)  , 
depend  on  the  initial  accelerometer  and  gravity  errors.  In 
the  alignment  process,  the  transformation  matrix  from  the 
platform  reference  frame  to  the  navigation  frame  is  rotated 
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into  alignment  with  the  sensed  gravity  vector.  This  causes 
the  initial  attitude  errors  to  correspond  to  errors  in  the 
sensed  gravity  vector.  The  alignment  heading  is  a  factor  in 
the  contribution  of  each  sensor  to  the  initial  state  errors 
(Ref  13) . 

Row  five  of  the  fundamental  matrix  can  be  used  to  show 
the  relationship  between  the  initial  east  tilt  error,  x  (7)  , 
and  acceleration  errors  (Ref  13) .  For  an  aircraft  at  rest, 
the  equation  is: 

0  .  -2  4?  sin  L  xq(4>  .  ( 7 )  xo<34>  CnJ[-  X0<3S)  Cny 

‘  Xo1431  Cny  fu  -  V491  <36» 

where 

Q  =  magnitude  of  the  earth's  angular  velocity 
f 

u*  force  up  in  the  body  frame 

Cnx  '  cny  '  cnz  *  elements  of  the  body  frame  to  navigation 

frame  transformation  matrix 

But  for  an  aircraft  at  rest,  fu=“fz=9  (where  fz  is  the  force 
in  the  z  direction  of  the  navigation  frame) ,  and  assuming 
alignment  at  a  random  heading  ,  the  equation  becomes: 

xq(7)  *  2£sin  L  xc(4)/g  +  cos  xQ(34)  -sin  (x0(35)/g 

-  xQ (43) )  -  xo(49)/g  (37) 

The  fourth  row  of  the  fundamental  matrix  can  be  used  to 
find  the  relationship  for  the  initial  north  tilt  error  xo(8) 


a>x=  ’Qh cos\l'  *  sin  i P 
(jjys  i2ucos  +  5U  sin i// 

Flight  Prof i le 

The  object  of  the  flight  profile  is  to  excite  the  long 
term  error  modes  of  the  system.  No  attempt  was  made  to 
follow  a  specific  mission  from  a  given  base  along  a  given 
flight  path,  but  rather  a  profile  was  developed  which 
included  representative  mission  segments  for  a  0130 
aircraft. 

Mission  Profile 

The  mission  profile  simulates  a  generic  C-130  mission. 
The  flight  includes  takeoff,  a  simple  departure,  a  high  level 
cruise  leg,  and  a  short  low  level  route.  The  20  segments  for 
this  mission  are  listed  in  Table  7  on  the  next  page.  In  the 
table,  time  is  given  in  seconds  duration  for  the  segment. 

The  acceleration  vector  is  divided  into  two  components;  one 
along  the  route  of  flight,  and  the  other  tangential  to  the 
flight  path  as  generated  by  the  aircraft  maneuvers. 

The  mission  starts  lined  up  on  the  runway  with  zero 
velocity.  The  start  point  is  35  degrees  north  latitude  and 
90  degrees  west  longitude  at  sea  level.  The  runway  heading 
is  an  arbitrary  315  degrees.  From  this  initial  condition, 
the  aircraft  accelerates  down  the  runway  until  it  reaches  a 
ground  speed  of  105  knots.  The  wind  is  assumed  to  be  zero. 

At  this  point,  the  aircraft  pitches  up  3  degrees  and  starts 
to  climb.  During  the  climb,  it  continues  to  accelerate  to  a 
climb  airspeed  of  180  knots.  The  departure  includes  two 
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heading  changes;  the  first  is  90  degrees  right  followed  by  a 
second,  at  level-off,  of  45  degrees  right.  The  entire 
departure  takes  just  under  33  minutes.  At  level-off,  the 
heading  is  straight  east  (heading  090  true) ,  and  the  aircraft 
begins  accelerating  to  a  cruise  airspeed  of  280  knots.  It 
should  be  noted  that  the  ground  speed  is  also  280  knots  since 
the  wind  is  assumed  to  be  calm.  The  high  altitude  cruise  leg 
is  exactly  6  hours  long. 

At  the  end  of  the  high  altitude  cruise  leg,  the  aircraft 

makes  a  90  degree  right  turn  to  a  heading  of  180  degrees  and 

starts  to  slow  to  its  low  altitude  cruise  airspeed  of  260 

knots.  At  the  descent  point,  the  aircraft  pitches  down  5.1 

degrees  and  descends  at  260  knots.  The  descent  takes  11 

minutes  and  continues  until  the  aircraft  levels  off  at  about 
900  feet.  The  low  level  portion  of  the  mission  has  10 

segments,  each  10  minutes  long.  The  low  level  segments  are 

separated  by  turns  of  45,  90,  or  135  degrees.  The  mission 

ends  at  low  level  after  a  total  flight  time  of  9  hours  2  1/2 

minutes . 

Summary 

Equation  1  expressed  the  form  of  a  set  of  linearized, 
stochastic,  first-order  differential  equations  which  can  be 
used  to  represent  an  error  model  of  the  SIGN-III  inertial 
navigation  system.  Equation  1  was  modified  by  the  addition 
of  two  different  models  of  the  vibration  environment  of  a 
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C-130A,  as  measured  during  flight  tests.  The  vibration 
models  represent  the  enviroment  through  a  series  of 
sinusoids.  While  an  accurate  representation  is  crucial,  of 
equal  importance  is  the  requirement  that  the  model  be 
computationally  efficient  since  the  vibration  to  be  modeled 
includes  frequencies  as  high  as  459  Hz.  In  addition,  flight 
mission  data  is  necessary  to  provide  inputs  to  the  system 
model . 
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Table  8.  Flight  Profile 


53 


III.  Software 


The  primary  program  .used  in  this  study  is  a  Monte  Carlo 
simulation  program  developed  by  the  Air  Force  Avionics 
Laboratory  called  SOFE  (Ref  10).  Two  additional  programs  used 
with  SOFE  are  a  flight  profile  generator,  PROFGEN  (Ref  9) ,  and 
a  statistical/plotting  postprocessor  for  SOFE,  SOFEPL  (Ref  4). 


SOFE 

This  section  will  discuss  the  implementation  of  SOFE  for 
the  system  simulation.  Appendix  A  has  a  discussion  of  the 
program  itself.  SOFE  can  be  used  to  implement  both  a  complete 
truth  model  and  a  reduced  order  Kalman  filter  model.  However, 
in  this  simulation  Kalman  filter  performance  was  not  an  issue 
and  the  filter  was  essentially  eliminated  by  making  its  order 
one,  the  minimum  allowable  in  SOFE. 

In  a  Monte  Carlo  simulation,  the  truth  model  is 
propagated  through  many  runs.  After  the  set  of  runs  is 
completed,  the  accumulated  data  can  be  analyzed  to  determine 
the  sample  mean  and  variance  of  the  states.  These  statistical 
computations  are  done  by  SOFEPL. 

SOFEPL 

SOFEPL  is  a  postprocessor  for  SOFE.  It  is  capable  of 
performing  statistical  computations  such  as  generating 
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ensemble  averages  and  standard  deviations  for  error  states. 

It  uses  a  graphics  package,  DISSPLA,  to  create  a  plot  file. 

The  DISSPLA  postprocessor  can  be  used  to  make  the  actual 
plots . 

PROFGEN 

Although  trajectory  information  can  be  generated  in 
SOFE,  a  separate  program,  PROFGEN,  was  used  to  simplify 
programming  and  to  reduce  the  amount  of  computer  resources 
reguired  for  any  given  run.  A  description  of  PROFGEN  is  given 
in  Appendix  B. 

PROFGEN  can  be  asked  to  generate  a  total  of  twenty 
output  variables,  but  only  seventeen  were  used  in  this 
simulation.  The  meanings  and  units  for  the  output  variables 
used  are  given  on  the  next  page.  Many  of  the  output  variables 
required  transformation  into  different  coordinate  frames 
before  they  could  be  used  in  SOFE.  PROFGEN  uses  a 
north-west-up  coordinate  frame  and  the  SOFE  simulation  was  run 
in  an  east-north-up  navigation  frame  and  a  fuselage 
bottom-nose-right  wing  body  frame.  The  transformation  from 
the  PROFGEN  frame  to  the  navigation  frame  was  easily  handled 
by  equating  the  corresponding  components  of  each  vector. 

Specific  force,  now  in  the  navigation  frame,  had  to  be 
transformed  into  the  body  frame  using  from  Equation  (2)  as 
shown: 
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Table  9.  PROFGEN  Output 


Symbol 

Definition 

Units 

t 

Time 

seconds 

L 

Latitude 

radians 

1 

Longitude 

radians 

h 

Altitude 

feet 

a 

Heading 

radians 

<t> 

Roll 

radians 

6 

Pitch 

radians 

V 

Yaw 

radians 

4> 

Roll  Rate 

rad/sec 

9 

Pitch  Rate 

rad/sec 

* 

Yaw  Rate 

rad/ sec 

vn 

North  Velocity  in  Nav  Frame 

ft/ sec 

-v 

West  Velocity  in  Nav  Frame 

ft/sec 

e 

vu 

Up  Velocity  in  Nav  Frame 

ft/ sec 

f  n 

Specific  Force  Along  North  Axis 

ft/sec 

-fe 

Specific  Force  Along  East  Axis 

ft/sec 

fu 

Specific  Force  Along  Up  Axis 

ft/sec 

C 

C 

xe 

xn 

XU 

C 

C 

ye 

yn 

yu 

C 

C 

ze 

zn 

ZU 

The  angular  velocity  of  the  body  frame  with  respect  to 
inertial  space  must  also  be  computed  from  the  PROFGEN  outputs. 
It  can  be  found  as  the  sum  of  the  angular  velocity  of  the 
earth  with  respect  to  inertial  space  plus  the  angular  velocity 
of  the  body  frame  with  respect  to  the  earth.  The  results  of 
these  computations  in  the  body  frame  are 


*  e 
4* 


"x 

= 

“z 

J 

where 


U)  *  -v  /R 
e  n 

cJn  “  Ve/R  +  cos  L 

=  (v  tan  L)/R  *  S3  sin  L 
wu  e 


R  =  20925640  ft 


S3  =  7.2921151  x  10 


rad/sec 


Summary 

This  study  used  three  software  programs:  SOFE,  SOFEPL, 
and  PROFGEN.  SOFE  is  a  Monte  Carlo  simulation  program  which 
was  used  to  propagate  the  SIGN-111  state  equations.  A  flight 
profile  generator,  PROFGEN,  was  used  to  provide  data  about  the 
aircraft  dynamics  to  SOFE.  SOFEPL  was  used  as  a 
statistical/plotting  postprocessor  for  the  SOFE  outputs. 
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IV.  Simulation  Results 


Program  Validation 

The  implementation  of  the  basic  nine-by-nine  error 
matrix  will  be  validated  by  comparing  its  response  to  various 
initial  conditions  to  those  obtained  by  Widnall  and  Grundy 
(Ref  3) .  They  generated  plots  showing  the  response  of  the 
unstable,  unaided  basic  nine-by-nine  system  model  to  initial 
errors  in  position,  velocity,  and  attitude.  The  plots  in 
Reference  3  include  the  response  to  initial  errors  in 
latitude  of  one  arc  minute,  an  altitude  error  of  ten  feet, 
north,  east  and  up  velocity  errors  of  one  foot  per  second, 
and  north,  east  and  up  attitude  errors  of  one  arc  minute. 
These  plots  were  duplicated,  but  only  a  limited  number  are 
included  for  sake  of  brevity.  Figures  10  and  11  show  the 
effect  of  an  initial  latitude  error  on  north  attitude  and 
latitude.  Figures  12  and  13  are  plots  of  the  longitude  and 
east  attitude  errors  resulting  from  initial  errors  in  north 
and  up  attitudes,  respectively.  The  errors  induced  in 
latitude  and  east  velocity  by  a  one  foot  per  second  error  in 
east  velocity  are  shown  in  Figures  14  and  15. 

A  small  scale  was  used  in  Reference  12,  and,  in  order  to 
allow  for  direct  comparison,  during  this  study.  Due  to  the 
small  scale,  it  is  difficult  to  be  certain  that  the  plots 
agree  completely.  However,  there  is  no  discernable 
difference  in  the  plots. 
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Chapter  2  included  the  development  of  the  basic  error 
model  and  two  different  methods  of  simulating  the  C-130A 
vibration  environment.  In  this  section,  the  results  of  Monte 
Carlo  simulation  using  those  models  will  be  presented.  The 
two  different  methods  of  modeling  the  vibration,  the 
sinusoidal  series  vibration  model  and  the  experimental 
integrated  effect  model,  were  simulated  in  separate 
15-minute,  20-run  Monte  Carlo  simulations.  During  these 
simulations,  the  only  aircraft  motion  was  that  simulated  by 
the  vibration  models.  The  results  of  these  runs  are  shown  in 
Figures  16  through  25.  Subsequently,  the  second  simplified 
method  was  simulated  for  an  8-hour  mission  with  a  20-run 
Monte  Carlo  simulation.  The  results  of  these  runs  are 
compared  with  a  20-run  simulation  of  the  50-state  basic  error 
model  without  vibration  in  Figures  26  through  37. 

Vibration- Induced  Errors 

While  20  runs  of  a  15-minute  mission  do  not  provide  a 
complete  picture  of  the  vibration-induced  errors,  computer 
limitations  prevented  the  use  of  a  longer  simulation.  The 
20-run  simulation  of  the  sinusoidal  series  model  required 
nearly  4  hours  of  cpu  time  on  a  CDC  Cyber  74  computer.  The 
experimental  integrated  effect  model  required  significantly 
less  time,  but  still  used  approximately  one  hour  of  cpu  time. 

Over  the  15-minute  runs,  the  vibration-induced  errors, 
in  both  cases,  were  very  small.  After  15  minutes,  the 
standard  deviations  of  the  latitude  error  for  the  first 
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Figure  21  .  Altitude  Error  for  Integrated  Effect  Method 
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method  were  approximately  0.3  micro-arcminutes  or  about  0.026 
inches.  The  standard  deviation  of  the  latitude  error  for  the 
second  method  was  even  smaller,  about  0.005  inches.  The 
error  standard  deviations ' from  the  second  simulation  were 
smaller,  by  a  factor  of  about  4.  The  longitude  errors  were 
larger;  from  the  first  simulation,  the  standard  deviation  of 
longitude  errors  after  15  minutes  was  about  20 
micro-arcminutes  (1.5  inches).  For  the  second  method,  the 
errors  were  again  smaller;  the  standard  deviations  were  about 
3.5  micro-arcminutes  (about  0.25  inches). 

The  plots  for  latitude,  longitude,  and  altitude 
demonstrate  another  significant  difference  in  the  results 
produced  by  the  two  methods.  The  plots  for  the  sinusoidal 
model  have  a  nonzero  mean,  while  those  for  the  integrated 
effect  method  appear  to  have  a  near  zero  mean.  These 
differences,  and  the  approximately  four-fold  difference  in 
the  magnitudes  of  the  standard  deviations,  is  probably  due  to 
a  flaw  in  the  assumption  that  the  state  transition  matrix  is 
time  invariant  for  the  sample  period  used. 

On  the  other  hand,  while  the  errors  are  smaller  than 
expected,  they  can  be  compared  with  the  results  of  a  similar 
study  (Ref  5)  published  by  Fontana  in  1972.  He  studied 
strapdown  inertial  navigation  for  the  European  space  vehicle. 
In  his  study,  he  examined  the  effects  of  sinusoidal 
vibrations  on  system  accuracy.  In  his  study,  he  set  the 
g-squared  gyro  sensitivities  at  0.01  degrees/hour/g-squared , 
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while  in  this  study  they  were  generated  randomly  with  a 
gaussian  distribution  and  a  standard  deviation  of  0.07 
degrees/hour/g-squared.  His  study  also  differed  in  that  the 
errors  were  determined  only  at  the  10  minute  point  in  the 
flight  (orbit  injection)  and  he  used  only  a  single  sinusoid 
with  results  computed  for  differing  vibration  amplitudes.  He 
also  assumed  a  worst  case  in  which  the  angle  between  the 
vibration  axes  and  the  gyro  spin  axes  was  45  degrees  (Ref  5) . 

Some  of  the  results  of  Fontana's  study  are  shown  in  Table 
9  on  the  following  page.  Vibrations  as  small  as  those 
normally  experienced  in  the  C-130A  were  not  examined,  but  some 
comparisions  can  still  be  made.  In  his  study,  the  error 
standard  deviations  caused  by  0.1  g  vibrations  were  exactly 
two  orders  of  magnitude  smaller  than  those  caused  by  1.0  g 
sinusoidal  vibrations,  indicating  that  the  standard  deviations 
of  the  induced  errors  are  proportionate  to  the  square  of  the 
vibration  level.  This  can  be  confirmed  by  examining  the 
g-squared  sensitive  gyro  drift  terms  in  the  S1GN-III  error 
model.  Equations  (44),  (45),  and  (46),  which  are  shown  below, 
are  the  g-squared  sensitive  gyro  drift  terms  as  extracted  from 
the  state  equations  for  states  7,  8,  and  9,  respectively. 

Under  the  assumptions  made  previously  about  the  nature  of  the 
vibration,  these  are  the  only  terms  through  which  vibration 
affects  system  accuracy. 


x7  (t) 


C  v  v  x,  <t)  +  C  v  v 
ex  n  u  19  ey  e  u 


(t) 
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Table  10.  RSS  Errors  Due  to  Sinusoidal  Vibration  (Ref  5) 


Vibration 

Velocity  Error 

Position  Error 

Amplitude 

(Ft/Sec) 

(Feet ) 

(G's) 

X  Y  2 

X  Y  Z 

0.1 

.0013 

.0017 

.0013 

.3542 

.  3739 

.1345 

1.0 

.1359 

.1708 

.1359 

35.40 

37.36 

13.43 

Table  11.  Comparison  of  Study  Results  (Ref  5) 


Position  Error 

(Feet) 

Fontana  Study 

Sinusoidal 

Integrated  Effect 

(Scaled) 

Method 

Method 

.2286 

.0996 

.0350 

-  C  V  v  X  ft)  ♦  OT 

ez  e  n  21 

x0(t)  *=  -  C  v  v  x  (t)  ♦  C  v  v  x  (t) 

8  nx  n  u  19  ny  e  u  20 

-  C  V  v  x„,  (t)  +  OT 

nz  e  n  21 

x-(t)  =  -  C  v  v  X  (t)  +  C  v  v  x  (t) 

9  zx  n  u  19  zy  e  u  20 

-  C  V  v  x,, (t)  *  OT 

zz  e  n  21 


where 

X.,  (t)  ,  x  (t)  ,  and  x  (t)  *  attitude  error  states,  east, 

7  o  9 

north,  and  up,  respectively 

x19(t)'  X20(t)'  and  X21(t)  =  ^-squared  gyro 

drift  coefficients 

C . •  =  element  of  the  transformation  matrix  from  local 
ID 

level  to  SIGN-III  coordinates 

ve ,  v  ,  and  v^  =  vibration  in  the  east,  north,  and  up 

directions,  respectively 

OT  =  other  non-g-squared  sensitive  terms 
From  these  equations,  it  is  evident  that  vibration  in  one 
axis  of  the  local  level  frame  is  multiplied  by  that  in  another 
axis  in  each  of  the  vibration  terms.  As  a  result,  vibration- 
induced  errors  should  be  proportional  to  the  product  of  the 
magnitudes  of  the  vibration.  If  the  vibration  levels  in  each 
axis  are  equal,  as  they  were  in  Reference  5,  then  the  errors 
would  be  proportional  to  the  square  of  the  vibration 
magnitude.  The  standard  deviations  of  the  errors  will  be 
affected  in  similar  fashion. 

Since  the  vibration  levels  and  the  gyro  drift  g-squared 
sensitivities  discussed  in  Reference  5  are  not  the  same  as 
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those  used  in  this  study,  the  resulting  error  standard 
deviations  must  be  scaled  before  comparisons  are  made.  This 
scaling  factor  was  determined  as  follows: 

S  *  (C/Cf)  x  (V/Vf)  (47) 

where 

S  *  the  scaling  factor 

C  *  the  standard  deviation  of  the  gyro  g-squared 
sensitive  drift  coefficients. 

Cf=  the  coefficients  used  in  Fontana's  study 
Vf e  the  sum  of  the  squares  of  the  vibration  levels  used 
in  Fontana's  study 
V  *  X  A  B  ♦XAC  +  X  B  C 

where 

A,  B,  C  are  the  vibration  levels  used  in  this  study. 


The  C/Cf  term  in  Equation  (47)  compensates  for  the  difference 
in  the  g-squared  gyro  drift  sensitivities.  The  V/Vf  term  is 
used  to  correct  for  the  differences  in  the  magnitudes  of  the 
sinusoidal  vibrations. 

Table  10  shows  the  RSS  position  errors  caused  by  0.1  g 
vibration  levels  from  Fontana's  study  scaled  down  by  a  factor 
of  2.8  as  stated  above.  The  RSS  errors  from  this  study  are 
also  included.  The  errors  resulting  from  this  study  are 
smaller  than  those  from  Fontana's  study,  but  are  of  the  same 
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order  of  magnitude  for  the  first  vibration  model  and  more 
than  one  order  of  magnitude  smaller  for  the  second  model. 

Much  of  the  differences  in  the  results  can  be  attributed  to 
Fontana's  assumption  of  worst  case  orientation  of  the 
vibration.  Thus,  the  sinusoidal  series  method  of  simulating 
vibration  compares  favorably  with  the  results  of  Fontana's 
study.  However,  it  is  unlikely  that  Fontana's  assumption  of 
worst  case  conditions  can  account  for  the  more  than  order  of 
magnitude  difference  in  the  results  between  his  study  and 
those  for  the  experimental  integrated  effect  method. 
Nonetheless,  the  computational  burden  associated  with  using 
the  sinusoidal  series  method  makes  its  use  impractical  unless 
the  mission  profile  is  extremely  short.  Based  on  the  CPU 
time  required  to  simulate  the  15-minute  missions,  it  woud 
have  required  exclusive  use  of  a  Cyber  74  computer  for  over 
10  days  to  complete  the  simulation  of  an  8-hour  mission.  As 
a  result,  it  could  not  be  used  and  the  full  8-hour  mission 
simulations  had  to  use  the  integrated  effect  method  in  spite 
of  its  demonstrated  lack  of  accuracy.  The  relative  efficiency 
of  the  this  method  enabled  the  full  8-hour  mission  to  be 
simulated  in  just  under  2  hours.  Thus,  this  method,  while 
not  accurate,  is  at  least  usable. 

Comparison  of  Errors  Due  to  Vibration  and  Total  System  Errors 


Figures  26  through  37  show  the  results  of  a  20-run, 
8-hour  simulaton  of  the  entire  system  model  without  vibration 
and  similar  runs  with  vibration  generated  by  the  experimental 
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integrated  effect  model  with  all  instrument  error  sources 
zeroed  out  with  the  exception  of  g-squared  effects.  The 
plots  of  the  vibration  runs  are  at  the  top  of  the  pages  with 
the  total  system  plots  of- the  same  error  state  at  the  bottom 
of  the  page.  While  Fontana's  study  indicates  that  the  second 
vibration  model  is  not  as  accurate  as  the  first  model, 
examination  of  the  15-minute  runs  of  these  methods  indicates 
that  the  second  method  has  some  limited  utility  if  the 
results  are  scaled  up  by  a  factor  of  approximately  five.  The 
resulting  loss  of  accuracy  is  more  than  compensated  by  the 
resulting  128-fold  reduction  in  cpu  time  required  which  makes 
the  second  method  usable  in  cases  in  which  the  first  method 
is  too  computationally  burdensome.  The  accuracy  of  the 
integrated  effect  method  is  too  poor  to  allow  its  use  for 
anything  other  than  rough  order  of  magnitude  estimates  of  the 
effects  of  the  vibration.  However,  it  can  give  an  indication 
of  the  order  of  magnitude  of  the  vibration-induced  errors  and 
can  help  in  determining  if  further  investigation  is 
warranted . 

Since  development  of  both  methods  was  based  on  the  same 
set  of  facts  and  assumptions,  with  one  exception,  it  is 
probable  that  the  additional  assumption  required  by  the 
integrated  effect  model  is  the  cause  of  the  loss  of  accuracy. 
This  assumption  was  that  the  state  transition  matrix  could  be 
treated  as  time-invariant  over  the  2  second  sampling  period. 
In  view  of  the  resulting  loss  of  accuracy,  it  is  likely  this 


assumption  was  not  valid. 

However,  examination  of  the  results  of  the  full  model 
simulations  and  of  the  vibration  simulations  reveals  that, 
even  if  the  vibration  is  Scaled  up  by  the  stated  factor  of 
five  (note:  the  plots  shown  do  not  include  this  scaling)  , 
the  effects  of  vibration  are  still  several  orders  of 
magnitude  smaller  than  the  errors  caused  by  other  factors. 
Thus,  if  either  of  the  vibration  models  accurately  represent 
the  C-130A  vibration  environment,  then  airframe  vibration  is 
not  a  major  source  of  inertial  navigation  system  errors. 
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V. 


Conclusions  and  Recommendations 


Conclusions 

This  study  examined  the  effects  of  airframe  vibration  on 
the  accuracy  of  the  SIGN-111  inertial  navigation  system,  in 
doing  so,  it  was  necessary  to  make  several  assumptions. 

First,  it  was  assumed  that  the  error  model  given  in  Widnall 
and  Grundy  is  accurate  and  complete  except  as  noted  in  Chapter 
2.  It  was  also  assumed  that  the  vibration  in  each  axis  was 
independent  in  phase  from  that  in  the  other  axes,  and  that  a 
sinusoidal  representation  of  the  vibration  was  accurate  enough 
to  cause  minimal  error  in  the  study  results. 

Two  vibration  models  were  developed  and  the  results  of 
Monte  Carlo  simulations  of  these  models  are  presented.  The 
results  of  the  sinusoidal  series  method  compare  favorably  with 
the  results  of  a  previous  study.  Unfortunately,  this  method 
is  prohibitively  burdensome  computationally  in  long 
simulations  involving  high  frequency  vibrations.  It  required 
about  1.5  CPU  seconds  per  mission  second  per  Monte  Carlo  run 
to  simulate  this  model.  The  integrated  effect  method  produced 
errors  which  were  smaller  by  a  factor  of  about  five,  but  was 
much  more  computationally  efficient.  It  ran  about  128  times 
faster  than  the  first  method.  The  results  of  this  method  might 
be  usable  if  the  required  degree  of  accuracy  is  very  low  and 
the  results  are  multiplied  by  a  factor  of  five  or  if  a 
correcting  coefficient  is  added  to  correct  the  error 
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magnitudes.  However,  at  best,  the  integrated  effect  method 
is  capable  of  giving  order  of  magnitude  estimates  of  the 
vibration-induced  errors.  As  such,  it  should  be  used  only  in 
determining  if  further  investigation  of  the  vibrational 
effects  is  required. 

Since  development  of  both  methods  was  based  on  the  same 
set  of  facts  and  assumptions,  with  one  exception,  it  is 
probable  that  the  additional  assumption  required  by  the 
integrated  effect  model  is  the  cause  of  the  loss  of  accuracy. 
This  assumption  was  that  the  state  transition  matrix  could  be 
treated  as  time-invariant  over  the  2  second  sampling  period. 
In  view  of  the  resulting  loss  of  accuracy,  it  is  likely  this 
assumption  was  not  valid. 

Regardless  of  which  vibration  model  is  used,  the 
resulting  navigation  errors  are  very  small,  ranging  from  4  to 
6  orders  of  magnitude  smaller  than  the  total  effects  of  other 
system  error  sources. 

Recommendations 

This  study  is  only  a  small  step  in  understanding  the 
effects  of  airframe  vibration  on  inertial  navigation  system 
errors  and  it  fell  short  of  its  goals.  However,  the 
computational  benefits  of  the  integrated  effect  method  may 
merit  additional  investigation  of  this  technique.  In 
particular,  the  sampling  period  should  be  examined  to  see  if 
a  shorter  sampling  period  will  result  in  improved  accuracy. 

On  the  other  hand,  since  shortening  the  sampling  period  will 
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also  increase  the  computational  burden  and  reduce  the  utility 
of  the  method,  it  is  recommended  that  sampling  periods 
between  0.1  and  0.5  seconds  be  used  as  a  starting  point. 
Furthermore,  the  factor  of  five  correction  used  with  the 
integrated  effect  vibration  model  was  obtained  by  comparing 
the  results  of  one  set  of  runs  produced  by  the  integrated 
effect  model  to  a  similar  set  produced  by  the  sinusoidal 
vibration  model.  No  analytical  basis  for  this  correction 
factor  was  found.  Thus,  the  correction  factor  used  in  this 
study  may,  or  may  not,  be  valid  for  a  different  inertial 
navigation  system  or  a  different  vibration  environment. 

Future  studies  should  examine  the  two  models  to  determine  an 
analytical  basis  for  the  correction  factor. 

In  addition,  studies  should  be  accomplished  to  determine 
whether  the  vibration  environment  is  well  represented  by  a 
sinusoidal  model.  It  should  also  be  determined  to  what 
degree  the  vibration  in  each  axis  is  correlated  with  that  in 
the  other  axes. 

In  the  system  model  used  in  this  study, 
vibration-induced  errors  entered  the  system  only  through  the 
g-squared  sensitive  gyro  drift  terms.  While  it  would  seem 
that  sinusoidal  vibration-induced  accelerometer  errors  should 
be  extremely  small,  they  should  be  studied  to  determine  if 
they  are  significant. 

This  study  was  limited  to  examining  only  one  strapdown 
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inertial  navigation  system.  While  the  results  should  be 
applicable  to  all  systems  of  comparable  type  and  quality, 
they  may  not  be  applicable  to  very  precise  systems  used  in 
long  range  ballistic  missiles  and  strategic  bombers  such  as 
the  B-52  and  B-1B.  This  study  was  also  limited  in  that  no 
attempt  was  made  to  examine  the  effect  of  vibration  on 
ring-laser  gyro  systems.  Indepth  studies  of  the  effects  of 
vibration  on  these  kinds  of  systems  should  be  accomplished. 
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Appendix  A 


SOFE  ;  A  Genera  1  ized  Digi  t-a  1  Simulation 
for  Optimal  Fi Iter  Eva  1 uat ion  (Ref  5) 

SOFE  is  an  efficient  general  purpose  program  which  was 
developed  for  the  design  and  evaluation  of  Kalman  filters  for 
use  in  integrated  systems.  (Ref  5)  Although  this  thesis  did 
not  involve  filter  design,  SOFE  was  used  to  perform  the 
Monte  Carlo  simulation  of  the  Inertial  Navigation  system 
error  state  equations  because  the  software  was  readily 
available,  provided  the  necessary  numerical  precision  and 
efficiency,  and  included  all  required  simulation 
capabilities . 

SOFE  is  divided  into  two  modules,  basic  SOFE  and 
user-written  SOFE.  Basic  SOFE  contains  31  routines  which 
perform  I/O,  problem  setup,  run  setup,  numerical  integration, 
measurement  update,  run  termination,  and  problem  termination. 
User-written  SOFE  consists  of  9  FORTRAN  subroutines  which 
specify  the  problem  to  be  simulated.  These  subroutines  allow 
for  filter  state  feedback,  computation  of  filter  matrices, 
computation  of  derivatives,  simulation  of  measurements, 
trajectory  data,  etc.  SOFE  was  designed  to  be  efficient  in 
the  use  of  memory  and  CPU  time.  This  was  accomplished  by 
dense  packing  of  arrays  and  vectors,  avoiding  the  use  of 
double  subscripts,  and  exploiting  the  symmetry  and/or  sparse 
properties  of  some  of  the  matrices.  (Ref  5) 
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The  SOFE  truth  model  is  an  implementation  of  the  error 
state  equations.  Only  4  of  the  available  9  user-written 
subroutines  were  used  to  implement  the  truth  model.  A  fifth 
subroutine,  CONVRT,  was  added  for  the  purpose  of  converting 
the  units  of  user  inputs  to  units  required  for  implementation 
of  the  truth  model.  The  other  subroutines  used  were  USRIN, 
SNOYS,  XSDOT ,  and  TRAJ.  USRIN  is  called  only  once  by  SOFE 
and  is  used  to  initialize  the  problem.  In  this  simulation, 
it  was  used  to  read  in  numerous  constants  and  the  standard 
deviations  of  the  initial  error  states.  It  was  also  used  to 
call  CONVRT.  SNOYS  was  used  to  simulate  the  dynamic  driving 
noises  by  injecting  white  noise  into  the  solution  of  the 
state  differential  equations  of  the  truth  model  states  at 
user-specified  intervals.  XSDOT  is  used  to  compute  the 
homogeneous  part  of  the  derivative  of  the  truth  state  vector 
xs.  TRAJ  was  used  only  to  set  some  constants  and  to  input 
trajectory  data  from  PROFGEN . 

Additional  data  are  entered  through  a  list  called 
PRDATA  in  CDC  NAMELIST  FORMAT.  PRDATA  includes  40  parameters 
which  remain  constant  throughout  the  simulation  and  are  used 
to  specify  the  user's  problem,  control  I/O,  and  regulate 
numerical  integration.  Listings  of  all  user-written 
subroutines  are  included  in  Appendices  C,  D,  and  E. 
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Appendix  B 


PROFGEN :  A  Computer  Program  for  Generating  Flight  Profiles 

PROFGEN  (Ref  6}  is  a  computer  program  which  computes 
flight  path  data  for  an  aircraft  flying  a  specified  route 
over  an  ellipsoidal  earth.  The  information  provided  includes 
position  (geographic  latitude,  longitude,  and  altitude) , 
velocity,  attitude,  and  attitude  rates  of  change.  Velocity 
is  computed  with  respect  to  the  earth  and  is  relative  to  a 
local  vertical  (navigation)  frame.  Acceleration  is  the  sum 
of  the  velocity  rates  of  change,  gravity,  and  Coriolis 
effects.  Attitude  is  expressed  in  terms  of  the  Euler  angles 
between  the  path  frame  and  the  navigation  frame:  roll, 
pitch,  and  yaw. 

PROFGEN  models  a  point  mass  and  does  not  model  the 
aerodynamics  of  the  aircraft.  As  a  result,  the  coordinate 
frame  of  the  flight  path  is  always  coincident  with  the  body 
coordinate  frame.  The  flight  profile  is  composed  of  up  to  50 
flight  segments,  each  segment  accomplishing  a  single  manuever 
from  the  set  of  five  possible  manuevers.  The  maneuvers 
available  are: 


Climb  or  Dive 
Coordinated  Turns 


90 


Sinusoidal  Heading  Changes 
Straight  Flights 
Rol  Is 

PROPGEN  allows  for  the  simulation  of  various  types  of 
aircraft  since  the  user  specifies  path  acceleration  rates, 
and  if  necessary,  centrifugal  acceleration  rates,  for  each 
segment.  The  final  value  for  each  variable  in  a  segment  is 
retained  and  used  as  the  initial  value  for  the  next  segment. 
As  a  result,  an  uninterrupted  time  history  exists  for  each 
variable . 

The  earth  model  used  in  PROFGEN  is  a  perfect  ellipsoid 
with  values  for  eccentricity ,  semimajor  axis  length,  spin 
velocity,  and  gravitational  constant  based  on  the  DOD 
Geodetic  System  1972.  Modelling  of  the  earth's  gravity 
includes  the  effects  of  latitude  and  altitude  changes  and  has 
both  radial  and  level  components  (Ref  6) .  While  this  model 
is  ' ot  overly  precise,  it  was  deemed  accurate  enough  for  use, 
without  revision,  in  this  study. 


APPENDIX  C 


SOFE  SUBROUTINES  FOR  50  STATE  MODEL  WITHOUT  VIBRATION 


SUBROUTINE  SNOYS { IRUN ,NF ,NS ,NXTJ ,XF , XS , XTRAJ) 
DIMENSION  XF (NF) ,XS (NS) , XTRAJ (NXTJ) 

COMMON/ SNOI S / SDWSO ( 50 ) , SDWS  (10) ,SDWFO,SDWF 
COMMON/DI ST/DALT , DGE , DGN , DGU 
COMMON /TRJCON / RE , G , OMEGA , RK 1 , RK  2 
V=SQRT ( XTRAJ { 8 ) **2  *  XTRAJ<9)**2) 

DT=T-TOLD 
SRDT-SQRT (DT) 

STDEVsSDWS (1) *SQRT (2*DT*V/DALT) 

CORNOS=GAUSS ( 0 . 0 ,STDEV) 

XS  (  6 ) *XS ( 6 )  +  RK2*CORNOS 

SDEV= (SDQS ( 2 )  ) *SRDT 

XS  (10) =XS (10)  ♦  GAUSS (0.0, SDEV ) 

SDEV® (SOWS (3) ) *SRDT 
XS(11)=XS(11)  +  GAUSS (0.0, SDEV) 

SDEV® (SDWS (4  )  ) *SRDT 

XS  ( 1 2 ) -XS ( 12 )  ♦  GAUSS (0.0, SDEV) 

SDEV® (SDWS (5)  ) *SRDT 

XS  (34) *XS  (34)  +  GAUSS (0.0, SDEV) 

SDEV- (SDWS (6) ) *SRDT 

XS  (35) =XS(35)  ♦  GAUSS (0.0, SDEV) 

SDEV® (SDWS (7) ) *SRDT 

XS  (  36) =  XS (36)  ♦  GAUSS (0.0, SDEV) 

XS  (46) =XS (46)  ♦  CORNOS 

SDEV= (SDWS (8 )  ) *SRDT 

XS  (48) =XS (48)  +  GAUSS(0.0,SDEV) 

SDEV® (SDWS (9) ) *SRDT 

XS  (49) ®XS(49)  +  GAUSS(0.0,SDEV) 

SDEV® (SDWS (10) ) *SRDT 

XS  (50) =XS (50)  ♦  GAUSS (0.0, SDEV ) 

RETURN 

ENTRY  SNOYSO 
TOLD-T 
RETURN 
END 

SUBROUTINE  USRIN 

COMMON/SNOIS/SDWSO(50) ,SDWS(10) , SDWFO,SDWF 
COMMON / D I ST / D ALT , DGE , DGN , DGU 
DIMENSION  DIST(4) 

NAMELIST/ IDIMEN/DIST, NWS, NWSO, SDWSO, SDWS, SDWFO,SDWF 
READ ( 5 , 1  DIMEN ) 
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WRITE ( 6 , IDIMEN) 

DALT«DIST ( 1 ) 

DGE=DIST(2) 

DGN*DIST ( 3 ) 

DGU=DIST ( 4 ) 

CALL  CONVRT  (NWSO ,NWS , SDWSO, SDWS ) 

RETURN 

END 

SUBROUTINE  XFDOT ( IRUN , T ,NF , NS , NXTJ , XF , XS , XTRAJ , NTR , PF , XDOT) 
DIMENSION  XF(NF) ,XS(NS) ,XTRAJ(NXTJ) ,  PF {NTR) ,XDOT(NF) 

XDOT ( 1 ) *0 . 0 
RETURN 

ENTRY  XFDOTO 
XF  ( 1 ) *0.0 
RETURN 
END 

SUBROUTINE  XSDOT ( IRUN ,T ,NF , NS , NXTJ , XF , XTRAJ , XDOT) 

DIMENSION  XF(NF) ,XS(NS) , XTRAJ (NXTJ) , XDOT (NS) 

COMMON/ SNOIS/ SDWSO (50) , SDWS (10) ,SDWFO,SDWF 

COMMON / TRJCON / RE , G , OMEGA , RK 1 , RK  2 

COMMON / TRA J 1  /  V , RLAT , RLON , ALT 

COMMON/ TRAJ2/VE, VN,VU , FE , FN ,FU ,WE,WN , WU 

COMMON / TRA J 3 /CEX , CEY ,CEZ , CNX ,  CNY ,CNZ,CUX,CUY,CUZ 

COMMON / TRA J  4 / S I NLAT , COSLAT , TANLAT 

COMMON/ TRAJ5 /OMEGAN , OMEGAU , RHOE , RHON , RHOZ , RKZ 

RLAT=XTRAJ  ( 1 ) 

RLON=XTRAJ (2) 

ALT=XTRAJ ( 4 ) 

ROLL® XT RAJ ( 5 ) 

PITCH=XTRAJ (6) 

YAW=XTRAJ ( 7 ) 

VE=XTRAJ ( 8 ) 

VN  =  XTRAJ ( 9 ) 

VU= XTRAJ ( 10 ) 

FE=XTRAJ (11) 

FN=XTRAJ (12) 

FU=XTRAJ (13) 

DROLL® XTRAJ (14) 

DP ITCH® XTRAJ ( 15 ) 

DYAW=XTRAJ ( 16 ) 

V=SQRT (VE*  *2  +VN*  *  2 ) 

SINLAT*SIN  (RLAT) 

COSLAT =COS (RLAT) 

TANLAT *S I NLAT /COSLAT 
S I NLON  =  S I N ( RLON ) 

COSLON=COS (RLON) 

OMEGAN® OMEGA* COSLAT 
OMEGAU =OMEGA* SI NLAT 
SX  =  SIN  (YAW) 

SY=SIN (ROLL) 
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SZ*SIN (PITCH) 

CX-COS(YAW) 

CY*COS  (ROLL) 

CZ=COS (PITCH) 

RHOE=-VN/RE 

RHON*VE/RE 

RHOU = VE  *T ANLAT /RE 

WE=RHOE 

WN=RHON  +  OMEGAN 
W(J=RHOU  +  OMEGAU 
RKZ“VU /RE 

F42-2.0* (OMEGAN*VN  ♦  OMEGAU#VU)  +  RHON*VN/COSLAT**2 
F  4  3  «  RHOU  *  RHOE  ♦  RHON*RKZ 
F 4 4  =  -RHOE * T ANLAT  -  RKZ 

F52““2 . 0*OMEGAN8VE  -  RHON*VE/COSLAT**2 

F53»RHON*RHOU  -  RHOE*RKZ 

F63«=2.0*G/RE  -  RHON**2  -  RHOE**2 

F92-WN  ♦  RHOU*TANLAT 

CXE=SX 

CXN=SZ*CY 

CXU=CZ*CY 

CYE=-CY*SX 

CYN»SZ*SY*SX  CZ*CX 
CYUbCZ*SY#SX  -  sz*cx 
CZE=-CY*CX 

CZN=“CZ*SX  ♦  SZ*SY*CX 

CZUeCZ*SY*CX  +  SY*SX 

CEX=-CXE 

CEY*CYE 

CEZ  sCZE 

CNX=CXN 

CNY  =CYN  . 

CZY  *CZN 
CUX*CXU 
CUY=CYl) 

CUZ=CZU 

FX=CXN*FN  ♦  CXE*FE  +  CXU*FU 

FY=CYN*FN  ♦  CYE*FE  ♦  CYU*FU 

FZ=CZN*FN  +  CZE*FE  +  CZU*FU 

WX=CXN*WN  +  CXE*WE  ♦  CXU*WU  ♦  DYAW 

WY«CYN*WN  ♦  CYE*WE  +  CYU*WU  +  DROLL 

WZ=CZN*WN  ♦  CZE*WE  +  CZU*WU  ♦  DPITCH 

WXP=0.0 

WXM  =  0 . 0 

WYP  =  0 . 0 

WYM- 0 . 0 

WZP=0.0 

WZM=0 . 0 

IF  (WX  .GE.  0.0  )  WXP*WX 
IF  (WX  .LE.  0.0  )  WXMSWX 
IF  (WY  .GE.  0.0  )  WYP=WX 
IF  (WY  .LE.  0.0  )  WYM=WX 
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IF  (WZ  .GE.  0.0  )  WZP*WZ 
IF  (WZ  .LE.  0.0  )  WZM=WZ 

XDOT ( 1) =  XS (2) *RHOU /COSLAT  -  XS ( 3 ) *RHON/ ( RE*COSLAT) 

1  ♦  XS (4) / (RE*COSLAT) 

XDOT ( 2 ) =XS ( 3 ) *RHOE/RE  +  XS(5)/RE 
XDOT ( 3 ) =  XS ( 6 ) 

XDOT  (4)*XS(2)*F42  XS(3)*F43  +  XS(4)*F44 

1  ♦  XS ( 5 ) * (WU  ♦  OMEGAU)  -  XS(6)*(WN  +  OMEGAN) 

1  -  XS ( 8 ) *FU  +  XS ( 9 ) *FN  ♦  XS(34)*CEX 

1  ♦  XS ( 3 5 ) *CEY  ♦  XS ( 36 ) *CEZ 

1  ♦  XS (37) *CEX*FX  ♦  XS (38) *CEY*FY 

1  +  XS (39) *CEZ*FZ  -  XS (40) *CEX*FZ 

1  ♦  XS (41) *CEX*FY  +  XS (42) *CEY*FZ 

1  -  XS (43) *CEY*FX  -  XS ( 44 ) *CEZ* FY 

1  4  XS (45) *CEZ*FX  4  XS  ( 48 ) 

XDOT (5)aXS(2)*F52  4  XS(3)*F53  -  XS(4)*2.0*WU  -  XS(5)*RKZ 
1  4  XS ( 6 ) *RHOE  4  XS ( 7 ) *FU  -  XS(9)*FE 

1  4  XS ( 3 4 ) *CNX  4  XS ( 3 5 ) *CNY  4  XS(36)*CNZ 

1  4  XS (37) *CNX*FX  4  XS (38) *CNY*FY 

1  4  XS (39) *CNZ*FZ  -  XS (40) *CNX*FZ 

1  4  XS ( 4 1 ) *CNX*FY  +  XS (42) *CNY*FZ  -  XS ( 4 3 ) *CNY*FX 

1  -  XS (44) *CNZ*FY  4  XS (45) *CNZ*FX  +  XS(49) 

XDOT ( 6 ) =  -  XS  (  2 ) *  2 . 0*OMEGAU*VE  4  XS(3)*(F63  -  RK1) 

1  4  XS ( 4 ) *2 . 0*WN  -  XS (5) *2.0*RHOE  -  XS(6)*RK2 

1  -  XS ( 7 ) *FN  4  XS ( 8 ) *FE  +  XS(34)*CUX 

1  4  XS ( 3 5 ) *CUY  4  XS ( 3 6 ) *CUZ  4  XS ( 3 7 ) *CUZ *FX 

1  4  XS ( 38 ) *CUY*FY  +  XS ( 39) *CUZ*FZ 

1  -  XS (40) *CUX*FZ  4  XS (41) *CUX*FY 

1  4  XS  (  4  2 ) *CUY*FZ  -  XS (43) *CUY*FX 

1  -  XS ( 44 ) CUZ*FY  4  XS(45)*CUZ*FX 

1  4  XS ( 46 ) * (RK1  -  RK2  +  V/DALT) 

1  4  XS (47) * (RK1*ALT  4  RK2*VU) 

1  4  XS ( 50) 

XDOT (7) =  -  XS ( 3 ) *RHCE/RE  -  XS(5)/RE  4  XS(8)*WU 
1  -  XS ( 9 ) *WN  4  XS  (10) *CEX  4  XS(11)*CEY 

1  4  XS  (12) *CEZ  ♦  XS (13) *CEX*FX 

1  4  XS (14) *CEX*FZ  4  XS (15) *CEY*FY 

1  4  XS (16) *CEY*FZ  4  XS (17) *CEZ*FZ 

1  -  XS (18) *CEZ*FY  -  XS (19) *CEX*FX*FY 

1  4  XS (20) *CEY*FX*FZ  -  XS (21) *CEZ*FX*FY 

1  4  XS  (22) *CEX*WXP  4  XS ( 2  3 ) *CEX*WXM 

1  4. XS (24) *CEY*WYP  ♦  XS ( 2 5 ) *CEY*WYM 

1  4  XS (26) *CEZ*WZP  +  XS (27) *CEZ*WZM 

1  4  XS (28) *CEX*WZ  -  XS (29) *CEX*WY 

1  -  XS (30) *CEY*WZ  4  XS (31) *CEY*WX 

1  4  XS (32) *CEZ*WY  -  XS (33) *CEZ*WX 

XDOT ( 8 ) =  -  XS ( 2 ) *OMEGAU  -XS ( 3 ) *RHON/ RE  4  XS(4)/RE 
1  -  XS(7)*WU  4  XS ( 9 ) *WE  4  XS ( 1 0 ) *CNX 

1  4  XS (11) *CNY  +  XS ( 1 2 ) *CNZ  4  XS ( 1 3 ) *CNX*FX 

1  4  XS (14) *CNX*FZ  4  XS (15) CNY* FY 

1  4  XS (16) *CNY*FZ  +  XS (17) *CNZ*FZ 

1  -  XS  (18) *CNZ*FY  -  XS (19) *CNX*FY*FZ 
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1  +  XS (20) *CNY*FX*FZ  -  XS ( 21) *CNZ*FX*FY 

1  +  XS (22) *CNX*WXP  ♦  XS ( 2 3 ) *CNX*WXM 

1  +  XS (26) *CNZ*WZP  +  XS (27) *CNZ*WZM 

1  +  XS  (28) *CNX*WZ  -  XS (29) *CNX*WY 

1  -  XS (30) *CNY*WZ  ♦  XS (31) *CNY*WX 

1  XS(32) *CNZ*WY  -  XS  (33)  *CNZ*WX 

XDOT ( 9 )  =  XS  (  2 ) *F92  -  XS ( 3 ) *RHOU / RE  ♦  XS  ( 4 ) *TANLAT/RE 
1  4  XS ( 7 ) *WN  -  XS ( 8) *WE  ♦  XS(10)*CUX  ♦  XS(11)*CUY 

1  +  XS ( 12) *CUZ  ♦  XS  (13) *CUX*FX  *  XS ( 14 ) *CUX*FZ 

1  +  XS (15) *CUY*FY  ♦  XS (16) *CUY*FZ 

1  +  XS ( 17 ) *CUZ*FZ  -  XS  ( 18 ) *CUZ*FY 

1  -  XS ( 1 9 ) *CUX*FY*FZ  ♦  XS ( 20 ) *CUY*FX*FZ 

1  -  XS (21) *CUZ*FX*FY  ♦  XS ( 2 2 ) *CUX*WXP 

1  +  XS ( 2 3 ) *CUX*WXM  ♦  XS (24) *CUY*WYP 

1  +  XS  (25)  *CUY*WYM  +  XS  ( 2 6 )  *Cl)Z *WZP 

1  4  XS (27) *CUZ*WZM  *  XS (28) *CUX*WZ 

1  -  XS ( 29 ) *CUX*WY  -  XS (30) *CUY*WZ 

1  4  XS (31) *CUY*WX  ♦  XS (32) *CUZ*WY 

1  -  XS  (  3  3 )  *C(JZ  *WX 

DO  10  1=10,45 

10  XDOT (I) =0.0 

XDOT (46)—  -  XS  (46) *V/DALT 
XDOT (47)  =  0.0 

XDOT (48)=  -  XS(48)*V/DGE 
XDOT (49)=  -  XS ( 49 ) *V/DGN 
XDOT ( 50 ) =  -  XS  (50) *V/DGU 
RETURN 

ENTRY  XSDOTO 
ALPHA=XTRAJ(3) 

DO  25  1=1,6 

25  XS (I ) =GAUSS(0.0,SDWSO(I) ) 

DO  45  1=10,50 

45  XS ( I ) =GAUSS(0.0,SDWSO(I) ) 

RLAT=XTRAJ (1) 

OMEGA=7. 2921151E-5 
OMEGAN=OMEGA*COS (RLAT) 

OMEGAU=OMEGA*S IN (RLAT) 

RK1  =  3E-2 
RK2  =  3E-4 
G«32. 0881576 
RE=20925640 . 0 
RANHED=GAUSS (0.0,1.8138) 

CHEAD=COS (RANHEAD) 

SHEAD=SIN (RHEAD) 

CNX= -CHEAD 
CUX=SHEAD 
SUY=CHEAD 
CNY-SHEAD 

WX=CNX*OMEGAN  +  CUX*OMEGAU 
WY=CUY*OMEGAU  4  CNY*OMEGAU 

XS(7) =XS (4) *2 . *OMEGAU/G  -  XS ( 6 ) *  2 . *OMEGAN/G 
1  4  XS ( 3  6 ) /G  -  XS  (  4  5 )  +  XS(48)/G 
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XS{8) *XS(5) *2.*0MEGAU/G  -  XS { 6 ) *2 . "OMEGAN/G 
1  ♦  XS ( 36) /G  -  XS  ( 45)  ♦  XS(48)/G 

XS (9) * (  -  XS (5) /RE  ♦  XS (8) "OMEGAU  +  XS(12) 

1  ♦  XS (31) *CNY*WX) /OMEGAU 

RETURN 
END 


I  C* 


SUBROUTINE  TRAJO ( IRUN , T , NF ,NS ,  NXTJ , X F , XS , XTRAJ ) 

DIMENSION  XF (NF) ,XS(NS) , XTRAJ (NXTJ) 

COMMON / TRJCON / RE , G , OMEGA , RK 1 , RK  2 
INTEGER  I TITLE (60) ,IPRSET(19) ,IRTSET(19) 

NAMELIST/ PRDATA/NSEGT , TSTART , VTO , ROLLO , PI TCHO , HEADO , ALF AO 
1  GLATO , TLONO , ALTO , I PRNT , I R I TE , I PLOT , ROLRAT , ROLTC , 

1  LLMECH ,LUNIT , RELERR , ABSERR , I PRSET , I RTSET 

C  READ  AND  ECHO  TITLE  AND  INPUT  DATA  FROM  PROFGEN  (TAPE  3) 

READ (3)  ITITLE, TOD AY, CLOCK 

READ ( 3 )  NSEGT , TSTART , VTO , ROLLO , P I TC HO , HEADO , ALFAO 
1  GLATO , TLONO , ALTO , I PRNT , I RI TE , I  PLOT , ROLRAT , 

1  ROLTC, LLMECH, LUNIT, RELERR, ABSERR, I PRSET , IRTSET 

DO  10  1*1,16 
10  READ (3)  DUM 

IF  (IRUN  .NE.  1)  RETURN 
WRITE ( 6 , PRDATA) 

WRITE (6,100) (ITITLE (I ) ,1*1,6) , TODAY, CLOCK 

OMEGA =70292115E-5 

G=32. 12698510 

RE*20925640 . 0 

RKl*3.0E-2 

RK2-3.0E-4 

100  FORMAT  ( / / 5 X , "THE  ABOVE  DATA  WAS  USED  TO  CREATE  THE  TRAJECTORY 
♦USING  'PROFGEN':" 

*  //10X, "TRAJECTORY  TITLE:  ",6A10, 

*  /10X, TRAJECTORY  RUN  DATE  AND  TIME:  " A10 , 5X , A10 ) 

RETURN 

END 


SUBROUTINE  CONVRT (NWSO , NWS , SDWSO , SDWS) 

DIMENSION  SDWSO (NWSO) , SOWS (NWS) 

PI = ABS (ACOS (-1.0) ) 

DEGREES  TO  RADIANS  CONVERSION 
CONVl*PI/ 180.0 
HOURS  TO  SECONDS  CONVERSION 
CONV2*3600 . 0 

ARCMINUTES  TO  RADIANS  CONVERTION 
CONV3= (1.0/60.0) "CONV1 
ARCSECONDS  TO  RADIANS  CONVERSION 
CONV4* (1.0/60.0) *CONV3 
G'S  TO  FEET  PER  SECOND  SQUARED  CONVERSION 
CONV5=32 . 2 

MICRO  G'S  TO  FEET  PER  SECOND  SQUARED  CONVERSION 
CONV6=32 . 2E-6 

DEGREES  PER  HOUR  TO  RADIANS  PER  SECOND  CONVERSION 


S;  vV- 
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CONV7-CONV1/CONV2 

C  PARTS  PER  MILLION  CONVERSION 

CCNV8-1.0E-6 

C  MICRORADIANS  TO  RADIANS  CONVERSION 

CONV9* 1 . 0E-6 

C  G  SQUARED  TO  (FEET  PER  SECOND  SQUARED)  CONVERSION 

CONV10-32.2**2 

C  SQUARE  ROOT  HOUR  TO  SQUARE  ROOT  SECOND  CONVERSION 

CONV11-60.0 

C  NAUTICAL  MILES  TO  FEET  CONVERSION 

CONV12-6076. 10333 
DO  15  1-1,2 

15  SDWSO(I) -SDWSO(I) *CONV3 
DO  20  1-10,12 

20  SDWSO(I) -SDWSO(I) *CONV7 
DO  40  1=13,18 

40  SDWSO ( I ) *SDWSO ( I ) *CONV7 /CONV5 
DO  50  1-19,21 

50  SDWSO ( I ) -SDWSO ( I ) *CONV7/CONV5**2 
DO  60  1-22,27 

60  SDCWSO ( I ) -SDWSO ( I ) *CONV8 
DO  70  1=28,33 

70  SDWSO ( I ) -SDWSO ( I ) *CONV4 
DO  80  1*34,36 

80  SDWSO ( I ) -SDWSO ( I ) *CONV6 
DO  90  1-37,39 

90  SDWSO ( I ) -SDWSO ( I ) *CONV8 
DO  100  1=40,45 

100  SDWSO ( I ) =SDWSO ( IO*CONV6 
DO  110  1*48,50 

110  SDWSO ( I ) *SDWSO ( I ) *CONV6 
DO  120  1*2,4 

120  SDWS(I) -SDWS (I) *CONVl/ (CONV2*CONVll) 

DO  130  1*5,7 

130  SDWS ( I ) -SDWS ( I ) *CONV6 /CONV1 1 
DO  140  1*8,10 

140  SDWS (I) -SDWS (I) *CONV6 
RETURN 
END 
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APPFND1 X  D 


SOFE  SUBROUTINES  FOR  SINUSOIDAL  VIBRATION  RUN 


SUBROUTINE  SNOYS ( IRUN,NF ,NS ,NXTJ , XF , XS ,XTRAJ ) 
DIMENSION  XF(NF) ,XS(NS) ,XTRAJ(NXTJ) 

DIMENSION  PSI (9) ,PHI(9) , THETA ( 9 ) ,AX(9) ,AY(9) ,AZ(9) 
DIMENSION  AXY  (9)  ,AXZ (9)  ,AYZ (9) 

COMMON/ SNOIS/ FU ,CEX,CEY,CEZ ,  CNX , CNY , CNZ ,CUX,CUY,CUZ 
COMMON/SNOIS2/ FXX , FYY ,  FZZ 

dt*t-told 
WT=3  20*T 
FX=FXX 
FY=FYY 
FZ=FZZ 
DO  10  1=1,9 

FX=FX  ^  COS (WT* I  «■  PSI (I) ) *AX{I) 

FY  =  FY  ♦  COS (WT* I  ♦  PHI ( I ) ) *AY ( I ) 

FZ*FZ  ♦  COS (WT*I  ♦  THETA (I) ) *AZ (I) 

FXFY=FX*FY 

FXFZ«FX*FZ 

FYFZ=FY*FZ 

VI*  -  DX*CEX*FYFZ  ♦  DY*CEY*FXFZ  -  D2*CEZ*FXFY 

V2  =  -  DX*CNX*FYFZ  ♦  DY*CNY*FXFZ  -  D2*CNZ*FXFY 

V3  =  -  DX*CUX*FYFZ  DY*CUY*FXFZ  -  DZ*CUZ*FXFY 

XS  ( 7 ) *XS ( 7 )  ♦  V1*DT 

XS ( 8 ) =XS ( 8 )  +  V2*DT 

XS (9) =XS (9)  ♦  V3*DT 

RETURN 

ENTRY  SNOYSO 

DO  30  1=1,9 

PHI ( I ) =GAUSS (0.0,1.81) 

PSI ( I ) *GAUSS (0.0,1.81) 

THETA ( I ) “GAUSS (0.0,1.81) 

30  CONTINUE 

AX(1) =0.2478 
AX(2) =0.2655 
AX ( 3 ) =0 . 2  2  3 
AX ( 4 ) =0 . 1 1 2 
AX(5) =0.204 
AX (6) =0 . 112 
AX(7)*0.129 
AX(8) =0.170 
AX(9) =0.091 
AY(1) =0.5657 
AY(2)=0.2325 
AY (3 ) =0 . 273 
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AY  (4) *0.170 
AY ( 5 ) -0.011 
AY  (6) -0.112 
AY  (7) -0.129 
AY ( 8 ) *  0 . 1 5  8 
AY ( 9 ) *  0 . 0  9 1 
A2 (1) *1.8238 
AZ  (2) *0.864 
AZ  ( 3 ) *0 . 8 15 
AZ  (4) *0. 288 
AZ  (5) *0 . 407 
AZ  (6) *0.500 
AZ  ( 7) *0 . 644 
AZ  (8) *0. 288 
AZ  (9) *0 . 204 
DO  40  1-1,9 

DX “GAUSS (0.0,3.273E-10) 

DY*GAUSS (0.0,3.273E-10) 

DZ “GAUSS (0.0,3. 273E-10) 

TOLD-T 

RETURN 

END 


SUBROUTINE  USRIN 

RETURN 

END 


SUBROUTINE  XFDOT ( I RUN ,T ,NF ,NS ,NXTJ , XF , XS , XTRAJ , NTR , PF , XDOT) 
DIMENSION  XF(NF) ,XS(NS) , XTRAJ (NXTJ) ,PF (NTR) ,XDOT(NF) 

XDOT(l) -0.0 
RETURN 

ENTRY  XFDOTO 
XF  (1) *0 . 0 
RETURN 
END 
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SUBROUTINE  XSDOT ( IRUN,T,NF,NS ,NXTJ ,XF , XTRAJ, XDOT) 
DIMENSION  XF(NF) ,XS(NS) , XTRAJ (NXTJ)  , XDOT  (NS) 

COMMON / T  RJCON / RE , G , OMEGA , RK 1 , R  K  2 
COMMON/ SNOIS2 /FXX , XYY, FZZ 

COMMON / SNOl S / FU ,CEX ,CEY , CEZ ,CNY ,CNY ,CNZ ,CUX ,CUY ,CUZ 
XDOT ( 1 ) “  XS (4) /  (RE*COSLAT) 

XDOT (2)*  XS ( 5 ) /RE 
XDOT (3) =  XS ( 6 ) 

XDOT ( 4 ) “  XS ( 5 ) *  2.0  *  OMEGAU  -  XS ( 6 ) *  2.0  *  OMEGAN 
1  -  XS (8) *  G 

XDOT (5) «  -  XS ( 4 ) *2 . 0*WU 

XDOT ( 6 ) “  XS(3) * (2.0*G/RE-RK1)  +  XS ( 4 ) *2 . 0*OMEGAN 
1  -  XS ( 6 ) *RK2 

XDOT (7) =  -  XS ( 5 ) /RE  ♦  XS(8)*OMEGAU  -  XS(9)*OMEGAN 
XDOT ( 8 ) *  -  XS (2) *OMEGAU*  ♦  XS(4)/RE  -  XS(7)*OMEGAU 
XPOT (9 ) =  XS ( 2 ) *OMEGAN  +  XS (4 ) *TANLAT/RE  +  XS(7)*OMEGAN 
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RETURN 

ENTRY  XSDOTO 
ROLL*XTRAJ (5) 

PITCH=XTRAJ (6) 

YAW*XTRAJ(7) 

FU=XTRAJ (13) 

SX*SIN (YAW) 

SY*SIN (ROLL) 

SZ*SIN (PITCH) 

CX«COS (YAW) 

CY*COS (ROLL) 

CZ-COS (PITCH) 

CEX-SX 
CEY*  -CY*SX 
CEZ=  -CY*CX 
CNXeSZ*CY 

CNY*SZ*SY*SX  ♦  CZ*CX 
CNZ*  -CZ*SX  *SZ*SY*CX 
CUX*CZ*CY 

CUY«CZ*SY*CX  ♦  SZ*CX 
CUZ«CZ*SY*CX  ♦  SZ*SX 
FXX=CUX*FU 
FYY*CUY*FU 
FZZ-CUX^FU 
ALPHA=XTRAJ (3) 

DO  25  1*1,9 
25  XS(I)-0.0 

RLAT=XTRAJ (1) 

COSLAT*COS (RLAT) 

TANLAT*TAN (RLAT) 

OMEGA® 7 .2921151 E- 5 
OMEGAN«=OMEGA*COSLAT 
OMEGAU=OMEGA*SINLAT 
RK1*  3  . OE-2 
RK2*3.0E-4 
G=32. 0881576 
RE* 209256 40 . 0 
SDEV* 1 . 0 

DX*GAUSS (0.0,3.273E-10) 

DY*GAUSS(0.0,3.273E-10) 

DZ*GAUSS(0.0,3.273E-10) 

RETURN 

END 

SUBROUTINE  TRAJO ( I RUN , T ,NF , NS ,  NXT  J , XF , XS , XTRAJ ) 

DIMENSION  XF(NF) ,XS(NS) , XTRAJ (NXTJ) 

COMMON /TRJCON/ RE, G , OMEGA , RK1 , RK2 
INTEGER  I TITLE (60) ,IPRSET(19) ,IRTSET<19) 

NAMELIST/ PRDATA/NSEGT , TSTART , VTO , ROLLO , P I TCHO , HE ADO , ALF AO 
1  GLATO , TLONO , ALTO , I PRNT , I RITE , I PLOT , ROLRAT , ROLTC , 

1  LLMECH , LUN IT , RELERR , ABSERR , I P "SET , I RTSET 

C  READ  AND  ECHO  TITLE  AND  INPUT  DATA  FROM  PROFGEN  (TAPE  3) 
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READ ( 3 )  ITITLE, TODAY, CLOCK 

READ ( 3 )  NSEGT , TSTART , VTO , ROLLO , P ITCHO , HEADO , ALFAO 
1  GLATO,TLONO, ALTO, IPRNT , I  RITE, I  PLOT , ROLRAT , 

1  ROLTC , LLMECH , LUN I T , RELERR, ABSERR, IPRSET , IRTSET 

DO  10  1*1,16 
10  READ (3)  DUM 

IF  (IRUN  .NE.  1)  RETURN 

UUTTP/fi 

WRITE  (6  \  100)  *(  ITITLE  (I)  ,1*1,6)  ,  TODAY, CLOCK 

OMEGA*70292115E-5 

G«32. 12698510 

RE-20925640.0 

RK1*3 . OE-2  N 

RK2*3 . OE-4 

100  FORMAT  (//5X, "THE  ABOVE  DATA  WAS  USED  TO  CREATE  THE  TRAJECTORY 
* USING  ' PROFGEN ' : " 

*  //10X, "TRAJECTORY  TITLE:  ”,6A10, 

*  /10X, TRAJECTORY  RUN  DATE  AND  TIME:  "A10,5X,A10) 

RETURN 

END 
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APPENDIX  E 


SOFE  SUBROUTINES  FOR  50  STATE  MODEL  WITH  VIBRATION 


SUBROUTINE  SNOYS ( IRUN ,NF ,NS ,NXTJ , XF , XS , XTRAJ) 
DIMENSION  XF(NF) ,XS(NS) , XTRAJ  (NXTJ) 

DIMENSION  PSI (9) ,PHI (9) , THETA (9) ,AX(9) ,AY(9) ,AZ (9) 
DIMENSION  AXY(9) ,AXZ (9) ,AYZ (9) 

COMMON / SNO I S / SDWSO (50) ,SDWS(10) ,SDWFO,SDWF 

COMMON/TRAJ2/VE,VN,VU,FE,FN,FU,WE,WN,WU 

COMMON /TRAJ 3 /CEX  ,  CEY ,  CEZ /CNX ,  CNY ,  CNZ ,CUX ,CUY ,CUZ 

COMMON / D I S  T / D ALT ,DGE,DGN,DGU 

COMMON / TRJCON / RE / G , OMEGA , RK 1 , RK2 

V*SQRT  ( XTRAJ  ( 8 )  *  *  2  <•  XTRAJ(9)**2) 

DT-T-TOLD 

SRDT*SQRT (DT) 

STDEV-SDWS (1) *SQRT(2*DT*V/DALT) 

CORNOS=GAUSS (0.0/ STDEV) 

XS  ( 6 ) *XS ( 6 )  ♦  RK2*CORNOS 
SDEV* (SDQS (2) ) *SRDT 
XS(10)-XS(10)  ♦  GAUSS (0.0, SDEV) 

SDEV* (SDWS (3) ) *SRDT 

fjf  XS(11)-XS(11)  ♦  GAUSS  (0.0,  SDEV) 

SDEV* (SDWS ( 4 ) ) *SRDT 

XS ( 12 ) “XS (12)  4  GAUSS (0.0, SDEV) 

SDEV* (SDWS (5) ) *SRDT 

XS  (34) *XS (34)  ♦  GAUSS (0.0, SDEV) 

SDEV* (SDWS (6) ) *SRDT 

XS (35) *XS ( 35)  +  GAUSS (0.0, SDEV) 

SDEV* (SDWS (7) ) *SRDT 

XS  (36) «XS (36)  +  GAUSS (0.0 ,SDEV) 

XS  (46) *XS (46)  ♦  CORNOS 

SDEV* (SDWS (8) ) *SRDT 

XS(48) *XS (48)  ♦  GAUSS (0.0, SDEV) 

SDEV* (SDWS (9) ) *SRDT 

XS (49) «XS (49)  +  GAUSS (0.0, SDEV) 

SDEV- (SDWS (10) ) *SRDT 

XS(50) *XS (50)  +  GAUSS (0.0, SDEV) 

DO  10  1*1,9 

PSI (I)-PSI (I)  ♦  GAUSS (0 . 0 , PHASE) 

PHI ( I ) *PHI (I)  ♦  GAUSS (0.0, PHASE) 

THETA (I) -THETA (I)  ♦  GAUSS (0 . 0 , PHASE) 

10  CONTINUE 
DX«XS  (19) 

DY«XS  (20) 

DZ-XS  (21) 
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DXCEX*DX*CEX 

DYCEY«DY*CEY 

DZCEZ«DZ*CEZ 

DXCNX*DX#CNX 

DYCNY=DY#CNY 

DZCNZ*DZ*CNZ 

DXCUX*DX*CUX 

DYCUY*DY*CUY 

DZCUZ«DZ*CUZ 

Vl-0.0 

V2-0 . 0 

V3-0 . 0 

DO  20  1*1/9 

V1*V1  +  (  - 


V3-V3  +  (  - 


V1*V1  +  (  -  DXCEX  * 

♦  DYCEY  * 
-  DZCEZ  * 

V2*V2  ♦  (  -  DXCNX  * 

♦  DYCNY  * 

.  -  DZCNZ  * 

V3-V3  +  (  -  DXCUX  * 

L  +  DYCUY  * 

L  -  DZCUZ  * 

CONTINUE 

XS ( 7 ) -XS ( 7 )  ♦  V1*DT 
XS  ( 8 )  *XS  ( 8 )  •»  V2*DT 
XS ( 9 ) *XS ( 9 )  4  V3*DT 
RETURN 

ENTRY  SNOYSO 

DO  30  1*1,9 

PHI (I) -GAUSS (0.0,1. 

PSI ( I ) *GAUSS (0.0,1. 

THETA (I) -GAUSS (0.0, 

CONTINUE 

AX(l)-0.2478 

AX(2)-0.2655 

AX ( 3 ) =0 . 223 

AX(4)*0.112 

AX(5)*0.204 

AX (6) *0 . 112 

AX ( 7} *0 . 129 

AX(8)-0.170 

AX  ( 9 ) *0 . 091 

AY(1)*0.5657 

AY(2)-0.2325 

AY ( 3) *0 . 273 

AY(4)*0.170 

AY  ( 5 )  ■  0 . 0 1 1 

AY  ( 6 )  *  0 . 1 1 2 

AY ( 7) -0 . 129 

AY ( 8 ) *0 . 158 

AY(9) *0.091 

AZ  (1) *1.8238 


AYZ  (I) 
AXZ  (I) 
AXY(I) 
AYZ  (I) 
AXZ  (I) 
AXY(I) 
AYZ (I) 
AXZ  (I) 
AXY(I) 


*  COS  (PHI (I) 

*  COS  (PSI (I) 

*  COS  (PSI (I) 

*  COS  (PHI ( I ) 

*  COS(PSI(I) 

*  COS  (PSI ( I ) 

*  COS (PHI (I) 

*  COS(PSI(I) 

*  COS  (PSI (I) 


-THETA ( I ) ) 

-  THETA (I) ) 
-PHI (I) ) ) / 2 
-THETA  (I) ) 

-  THETA (I) ) 
-PHI ( I ) )  )  / 2 
-THETA  ( I ) ) 

-  THETA (I) ) 
-PHI  (I) ) )/2 


81) 

81) 

1.81) 


AZ  ( 2 ) =0 . 864 
AZ  (3) =0.815 
AZ  ( 4 ) *0 . 28 8 
AZ  ( 5) *0 . 407 
AZ  (6) =0.500 
AZ  (7) =0.644 
AZ  (8) =0.288 
AZ  ( 9 ) =0 . 204 
DO  40  1=1,9 
AXY (I)«AX (I) *AY(I) 

AXZ ( I ) * AX ( I ) * AZ { I ) 

AYZ ( 1 ) ■ AY { I ) * AZ ( I ) 

40  CONTINUE 
TOLD=T 
RETURN 
END 

SUBROUTINE  USRIN 

COMMON/ SNOIS/SDWSO (50) , SDWS (10) , SDWFO, SDWF 
COMMON / D I ST / D ALT , DGE , DGN , DG  U 
DIMENSION  DIST  ( 4 ) 

NAMELIST/ I DIMEN/DI ST , NWS , NWSO , SDWSO , SDWS , SDWFO , SDWF 
READ (5 , I DIMEN) 

WRITE (6 , I DIMEN) 

DALT  =  DIST  ( 1 ) 

DGE=DIST ( 2 ) 

DGN=DIST ( 3 ) 

DGU=DIST ( 4 ) 

CALL  CONVRT  (NWSO , NWS , SDWSO, SDWS ) 

RETURN 

END 

SUBROUTINE  XFDOT ( I RUN , T , NF , NS ,NXT  J , XF , XS , XTRA J , NTR , PF , XDOT ) 
DIMENSION  XF (NF ) ,XS(NS) , XTRAJ (NXTJ) , PF (NTR) ,XDOT(NF) 

XDOT (1) =0.0 
RETURN 

ENTRY  XFDOTO 
XF ( 1 ) =0.0 
RETURN 
END 

SUBROUTINE  XSDOT ( IRUN, T ,NF ,NS ,NXT J , XF , XTRAJ , XDOT) 

DIMENSION  XF (NF) ,XS(NS) , XTRAJ (NXTJ) , XDOT (NS) 

COMMON/ SNOIS/ SDWSO (50) , SDWS (10) , SDWFO, SDWF 

COMMON /TRJCON / RE , G , OMEGA , RK 1 , R  K  2 

COMMON / TRA J 1  /  V , RL AT , RLON , ALT 

COMMON/ TRAJ 2 /VE ,VN ,  VU , FE , FN , FU , WE ,  WN , WU 

COMMON/ TRAJ 3  / CEX ,  CEY  ,CEZ  ,CNX  ,  CNY  ,  CNZ  ,CU,: ,  CUY ,  CUZ 

COMMON/TRAJ4/SINLAT, COSLAT, TANLAT 

COMMON /TRAJ  5 / OMEG AN , OMEGAU , RHOE , RHON , RHOZ , RK  Z 

RLAT-XTRAJ(l) 

RLON=XTRAJ ( 2 ) 
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ALT*XTRAJ (4 ) 

R0LL«XTRAJ(5) 

PITCH*XTRAJ (6) 

YAW*XTRAJ ( 7 ) 

VE=XTRAJ ( 8 ) 

VN=XTRAJ { 9 ) 

VU*XTRAJ ( 10 ) 

FE=XTRAJ (11) 

FN=XTRAJ (12) 

FU*XTRAJ (13) 

DROLL*  XTRA  J (14) 

DPITCH*XTRAJ ( 15 ) 

DYAW*XTRAJ (16) 

V*SQRT (VE**2+VN**2) 

SINLAT«SIN{RLAT) 

COSLAT*COS ( RLAT ) 

TANLAT*SINLAT /COSLAT 
SINLON=SIN ( RLON ) 

COSLON*COS ( RLON ) 

OMEGAN«OMEGA*COSLAT 
OMEGAU “OMEGA* S I NLAT 
SX“SIN (YAW) 

SY“SIN (ROLL) 

SZ-SIN (PITCH) 

CX*COS ( YAW) 

CY  *COS (ROLL) 

CZ“COS (PITCH) 

RHOE*-VN/RE 

RHON=VE/RE 

RHOU  = VE*  T ANLAT / RE 

WE«=RHOE 

WN=RHON  ♦  OMEGAN 
WU=RHOU  ♦  OMEGAU 
RKZ  * VU / RE 

F42*2 . 0* ( OMEGAN *VN  ♦  OMEGAU*VU)  ♦  RHON*VN/COSLAT* * 2 
F43“RHOU*RHOE  ♦  RHON*RKZ 
F44=-RHOE*TANLAT  -  RKZ 

F52*-2 . 0*OMEGAN8VE  -  RHON*VE/COSLAT* *  2 

F53=RHON*RHOU  -  RHOE*RKZ 

F63*2 . 0*G/RE  -  RHON**2  -  RHOE**2 

F92“WN  *  RHOU*TANLAT 

CXE*SX 

CXN“SZ*CY 

CXU=CZ*CY 

CYE“-CY*SX 

CYN*SZ*SY*SX  ♦  CZ*CX 
CYU=CZ*SY*SX  -  SZ*CX 
CZE“-CY*CX 

CZN“-CZ*SX  ♦  SZ*SY*CX 
CZU“CZ*SY*CX  +  SY*SX 
CEX=CXE 
CEY*CYE 
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CEZ«CZE 

CNX-CXN 

CNY=CYN 

CZY-CZN 

C(JX*CXU 

CUY-CYU 

CUZ*CZU 

FX-CXN*FN  ♦  CXE*FE  ♦  CXU*FU 

FY*CYN*FN  ♦  CYE*FE  ♦  CYU*FU 

FZ«CZN*FN  +  CZE*FE  ♦  CZU*FU 

WX*CXN»WN  ♦  CXE*WE  ♦  CXU*WU  ♦  DYAW 

WY*CYN*WN  +  CYE*WE  +  CYU*WU  +  DROLL 

WZ«CZN*WN  ♦  CZE*WE  ♦  CZU*WU  +  DPITCH 

WXP'0.0 

WXM=0 . 0 

WYP<=0.0 

WYM*0 . 0 

WZP*0.0 

WZM=0 . 0 

IF  (WX  .GE.  0.0  )  WXP«WX 
IF  (WX  .LE.  0.0  )  WXM*WX 
IF  (WY  .GE.  0.0  )  WYP=WX 
IF  (WY  .LE.  0.0  )  WYM*WX 
IF  (WZ  .GE.  0.0  )  WZP=WZ 
IF  (WZ  .LE.  0.0  )  WZM=WZ 

XDOT (1) -XS (2) *RHOU/COSLAT  -  XS ( 3 ) *RHON/ { RE*COSLAT) 

L  +  XS ( 4 ) / (RE*COSLAT ) 

XDOT ( 2 ) =XS ( 3 ) *RHOE/RE  ♦  XS(5)/RE 
XDOT ( 3 ) -XS (6) 

XDOT ( 4 ) KXS ( 2 ) *F4  2  ♦  XS(3)*F43  +  XS(4)*F44 
1  ♦  XS(5)*(WU  ♦  OMEGAU)  -  XS(6)*(WN  ♦  OMEGAN) 

1  -  XS (8) *FU  ♦  XS ( 9 ) *FN 

1  ♦  XS(34)*CEX  ♦  XS ( 35} *CEY  ♦  XS(36)*CEZ 

1  ♦  XS (37) *CEX*FX  ♦  XS (38) *CEY*FY 

1  ♦  XS ( 39) *CEZ*FZ  -  XS  ( 40 ) *CEX*FZ 

1  +  XS (41) *CEX*FY  +  XS { 42 ) *CEY*FZ 

1  -  XS ( 4 3 ) *CEY*FX  -  XS ( 44 ) *CEZ*FY 

1  ♦  XS (45) *CEZ*FX  +  XS ( 4  8 ) 

XDOT ( 5) =XS ( 2 ) *F52  ♦  XS(3)*F53  -  XS(4)*2.0*WU  -  XS(5)*RKZ 
1  ♦  XS ( 6 ) *RHOE  ♦  XS  (  7 ) *FU  -  XS(9)*FE 

1  +  XS ( 34) *CNX  +  XS(35)*CNY  ♦  XS(36)*CNZ 

1  ♦  XS (37) *CNX*FX  ♦  XS ( 38 ) *CNY*FY 

1  ♦  XS (39) *CNZ*FZ  -  XS ( 4 0 ) *CNX*FZ 

1  +  XS (41) *CNX*FY  ♦  XS (42) *CNY*FZ 

1  -  XS (43) *CNY*FX  -  XS ( 4 4 ) *CNZ*FY 

1  ♦  XS ( 45) *CNZ*FX  ♦  XS (49) 

XDOT ( 6 ) ■  -  XS ( 2 ) *  2 . 0*OMEGAU*VE  ♦  XS(3)*(F63  -  RK1 ) 

1  ♦  XS ( 4) *2 . 0*WN  -  XS ( 5 ) *2 . 0*RHOE 


1  -  XS ( 6 ) *RK2  -  XS ( 7 ) *FN 

1  +  XS ( 8 ) *FE  ♦  XS ( 34 ) *CUX  ♦  XS(35)*CUY 

1  4  XS ( 36) *CUZ  *  XS ( 37) *CUZ*FX  +  XS ( 3 8 ) *CUY*FY 

1  4  XS (39) *CUZ*FZ  -  XS (40) *CUX*FZ  4  XS (41) *CUX*FY 


+ 

+ 


XS (22) *CEX*WXP 
XS (24) *CEY*WYP 
XS  (26) *CEZ*WZP 
XS  (28) *CEX*WZ 
XS (30) *CEY*WZ 
XS (32) *CEZ*WY 


1  ♦  XS ( 42 ) *CUY*FZ  -  XS (43) *CUY*FX 

1  -  XS (44 ) CUZ*FY  ♦  XS (45) *CUZ*FX 

1  XS  (46 )  *  (RK1  -  RK2  ♦  V/DALT) 

1  ♦  XS (47) * (RK1*ALT  ♦  RK2*VU) 

1  +  XS ( 50 ) 

XDOT ( 7 ) *  -  XS ( 3 ) *RHOE/RE  -  XS(5)/RE  +  XS(8)*WU 
XS ( 9 ) *WN  ♦  XS ( 10 ) *C.EX  ♦  XS(11)*CEY 
XS  (12) *CEZ  +  XS ( 13 ) *CEX*FX 
XS (14 ) *CEX*FZ  ♦  XS ( 15 ) *CEY*FY 
XS  (16) *CEY*FZ  ♦  XS (17) *CEZ*FZ 
XS (18) *CEZ*FY  -  XS (19) *CEX*FX*FY 
XS  (20) *CEY*FX*FZ  -  XS  (21) *CEZ*FX*FY 

XS (23) *CEX*WXM 
XS (25) *CEY*WYM 
XS (27) *CEZ  *WZM 
XS (29) *CEX*WY 
♦  XS (31) *CEY*WX 
-  XS (33) *CEZ*WX 

XDOT (8) =  -  XS ( 2 ) *OMEGAU  -XS ( 3 ) *RHON/RE  4  XS(4)/RE 

-  XS ( 7 ) *WU  *  XS ( 9 ) *WE  +  XS ( 1 0 ) *CNX 
4  XS ( 1 1 ) *CNY  4  XS ( 1 2 ) *CNZ 
4  XS (13) *CNX*FX  +  XS  (14) *CNX*F2 

4  XS ( 1 5 ) CNY*FY  4  XS (16) *CNY*FZ 
4  XS ( 1 7 ) *CNZ *FZ  -  XS  (18) *CNZ*FY 

-  XS (19) *CNX*FY*FZ  4  XS  (  20) *CNY*FX*FZ 

-  XS ( 2 1 ) *CNZ*FX*FY  ♦  XS (22) *CNX*WXP 
4  XS (23) *CNX*WXM  4  XS (26) *CNZ*WZP 
4  XS (27) *CNZ*WZM  +  XS (28) *CNX*WZ 

-  XS (29) *CNX*WY  -  XS (30) *CNY*WZ 
♦  XS ( 3 1 ) *CNY*WX  +  XS ( 3  2 ) *CNZ*WY 

-  XS (33) *CNZ*WX 

XDOT ( 9 )  =  XS ( 2 ) *F92  -  XS ( 3 ) *RHOU/RE  4  XS ( 4 ) *TANLAT/RE 
1  ♦  XS ( 7 ) *WN  -  XS ( 8 ) *WE  ♦  XS(10)*CUX  4  XS(11)*CUY 

1  4  XS  ( 1  2  )  *CUZ  4  XS  (13)  *CUX*FX  ♦  XS  ( 1 4 ) '»CUX*FZ 

1  4  XS (15) *CUY*FY  4  XS (16) *CUY*FZ 

1  4  XS(17) *CUZ*FZ  -  XS(18)*CUZ*FY 

1  -  XS ( 1 9 ) *CUX*FY*FZ  +  XS  (20) *CUY*FX*FZ 

1  -  XS (21) *CUZ*FX*FY  +  XS (22) *CUX*WXP 

1  4  XS ( 2 3 ) *CUX*WXM  4  XS (24) *CUY*WYP 

1  4  XS ( 25 ) *CUY*WYM  ♦  XS { 26 ) *CUZ*WZP 

1  4  XS  (  27  )  *Cl)Z*WZM  4  XS  (28)  *CUX*WZ 

1  -  XS ( 29 ) *CUX*WY  -  XS (30) *CUY*WZ 

1  4  XS (31) *CUY*WX  4  XS (32) *CUZ*WY 

1  -  XS ( 3  3 ) *CUZ*WX 


1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 


DO  10  1=10,45 
10  XDOT (I) =0.0 

XDOT  (46)=  -  XS  (46) *V/DALT 
XDOT(47) =0.0 

XDOT (48)=  -  XS  (  4  8 ) *V/DGE 
XDOT ( 49  )  =  -  XS  (  4  9 ) *V/DGN 
XDOT ( 50 )  =  -  XS  (  50 ) *V/DGU 
RETURN 
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ENTRY  XSDOTO 
ALPHA-XTRAJ (3) 

DO  25  1-1,6 

25  XS(I) -GAUSS(0.0,SDWSO(I)  ) 

DO  45  1-10,50 

45  XS(I) «GAUSS(0.0,SDWSO(I)  ) 

RLAT-XTRAJ (1) 

OMEGA-7. 2921151E-5 
OMEGAN -OMEGA *COS  (RLAT) 

OMEG AU • OMEG A  *  S I N ( RLAT ) 

RK1-3E-2 

RK2-3E-4 

G-32. 0881576 

RE-20925640.0 

RANHED-GAUSS (0.0,1.8138) 

CHEAD-COS (RANHEAD) 

SHEAD-SIN (RHEAD) 

CNX— CHEAD 
CUX-SHEAD 
SUY-CHEAD 
CNY-SHEAD 

WX«CNX*OMEGAN  +  CUX*OMEGAU 
WY=CUY*OMEGAU  +  CNY*OMEGAU 

XS(7)«XS(4)*2.*OMEGAU/G  -  XS ( 6 ) *2 . *OMEGAN/G 
1  +  XS (36) /G  -  XS (45 )  ♦  XS(48)/G 

XS ( 8 ) -XS ( 5 ) *2 . *OMEGAU/G  -  XS ( 6 ) * 2 . *OMEGAN/G 
1  ♦  XS ( 36 ) /G  -  XS (45)  ♦  XS(48)/G 

XS ( 9 ) - (  -  XS ( 5 ) /RE  ♦  XS ( 8 ) *GMEGAU  ♦  XS  { 1 2 ) 

1  ♦  XS ( 31 ) *CNY*WX) / OMEGAU 

RETURN 
END 

SUBROUTINE  TRAJO ( I RUN ,T,NF,NS, NXTJ , XF , XS , XTRAJ) 

DIMENSION  XF (NF)  ,XS (NS) , XTRAJ (NXTJ) 

COMMON / TRJCON / RE , G , OMEGA , RK 1 , RK  2 
INTEGER  ITITLE (60) ,IPRSET(19) ,IRTSET(19) 

NAMELIST/PRDATA/NSEGT,TSTART,VTO,ROLLO,PITCHO,HEADO,ALFAO 
1  GLATO , TLONO , ALTO ,IPRNT,IRITE,I PLOT , ROLR AT , ROLTC , 

1  LLMECH,LUNIT,RELERR, ABSERR, IPRSET, IRTSET 

C  READ  AND  ECHO  TITLE  AND  INPUT  DATA  FROM  PROFGEN  (TAPE  3) 
READ (3)  ITITLE, TODAY, CLOCK 

READ ( 3 )  NSEGT , TSTART , VTO , ROLLO , P I TCHO , HEADO , ALFAO 
1  GLATO , TLONO , ALTO ,  I PRNT , I R I TE , I PLOT , ROLRAT , 

1  ROLTC, LLMECH,LUNIT,RELERR, ABSERR, IPRSE1 , IRTSET 

DO  10  1-1,16 
10  READ (3)  DUM 

IF  (IRUN  .NE.  1)  RETURN 
WRITE (6 ,PRDATA) 

WRITE (6, 100) (ITITLE (I) ,1-1,6) , TODAY, CLOCK 
OMEGA-70  2 92115E-5 
G-32. 12698510 
RE-20925640.0 
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RK1-3.0E-2 

RK2*3.0E-4 

100  FORMAT  ( / /5X , "THE  ABOVE  DATA  WAS  USED  TO  CREATE  THE  TRAJECTORY 

♦USING  ' PROFGEN ' : " 

*  //10X, "TRAJECTORY  TITLE:  ",6A10, 

♦  /10X, TRAJECTORY  RUN  DATE  AND  TIME:  "A10,SX,A10) 

RETURN 

END 

SUBROUTINE  CONVRT (NWSO,NWS , SDWSO , SDWS) 

DIMENSION  SDWSO (NWSO) , SDWS (NWS) 

PI*  ABS ( ACOS (-1.0) ) 

C  DEGREES  TO  RADIANS  CONVERSION 

C0NV1*PI/ 180.0 

C  HOURS  TO  SECONDS  CONVERSION 

CONV2-3600.0 

C  ARCMINUTES  TO  RADIANS  CONVERTION 

CONV3* (1.0/60.0) *CONVl 

C  ARCSECONDS  TO  RADIANS  CONVERSION 

CONV4* (1.0/60.0) *CONV3 

C  G'S  TO  FEET  PER  SECOND  SQUARED  CONVERSION 
CONV5*32 . 2 

C  MICRO  G'S  TO  FEET  PER  SECOND  SQUARED  CONVERSION 
CONV6-32.2E-6 

C  DEGREES  PER  HOUR  TO  RADIANS  PER  SECOND  CONVERSION 

CONV 7 *CONV 1 / CONV 2 

C  PARTS  PER  MILLION  CONVERSION 

CONV3* 1 . OE-6 

C  MICRORADIANS  TO  RADIANS  CONVERSION 

CONV9*l . OE-6 

C  G  SQUARED  TO  (FEET  PER  SECOND  SQUARED)  CONVERSION 
CONV10=32.2**2 

C  SQUARE  ROOT  HOUR  TO  SQUARE  ROOT  SECOND  CONVERSION 

CONV11=60.0 

C  NAUTICAL  MILES  TO  FEET  CONVERSION 

CONV12*6076. 10333 
DO  15  1*1,2 

15  SDWSO ( I ) *SDWSO ( I ) *CONV3 
DO  20  1*10,12 

20  SDWSO ( I ) *SDWSO ( I ) *CONV7 
DO  40  1*13,18 

40  SDWSO ( I ) =SDWSO ( I ) *CONV7 / CONV5 
DO  50  1*19,21 

50  SDWSO ( I ) *SDWSO ( I ) *CONV7/CONV5  **  2 
DO  60  1*22,27 

60  SDCWSO ( I ) *SDWSO ( I ) *CONV8 
DO  70  1*28,33 

70  SDWSO ( I ) “SDWSO ( I ) *CONV4 
DO  80  1*34,36 

80  SDWSO ( I ) *SDWSO ( I ) *CONV6 
DO  90  1*37,39 

90  SDWSO ( I ) *SDWSO ( I ) *CONV8 
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DO  100  1-40,45 

100  SDWSO ( I ) -SDWSO ( IO*CONV6 
DO  110  1-48,50 

110  SDWSO ( I ) -SDWSO ( I ) *CONV6 
DO  120  1-2,4 

120  SDWS (I) -SDWS (I) *CONVl/ (CONV2*CONVll) 
DO  130  1-5,7 

130  SDWS ( I ) -SDWS  < I ) *CONV6/CONVll 
DO  140  1-8,10 

140  SDWS ( I ) -SDWS ( I ) *CONV6 
RETURN 
END 
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